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ABSTRACT 


The  present  study  was  undertaken  in  an  effort  to  im¬ 
prove  numerical  models  for  meso-scale  and  small-scale  effects 
which  influence  global  weather  and  its  modification.  Two 
major  areas  are  being  studied:  the  effects  of  mountain 
ranges  on  energy  and  momentum  transfer,  and  the  transient 
interaction  of  solar  radiation  with  the  earth's  atmosphere. 

It  is  hoped  that  the  results  of  these  studies  will  lead  to 
calculationally  inexpensive  prescriptions  which  can  be  in¬ 
corporated  into  meso-scale  and  global-scale  atmospheric  cir¬ 
culation  codes. 


3SR-795 


TABLE  OF  CONTENTS 

Page 

ABSTRACT .  ii 

NOMENCLATURE  FOR  SECTIONS  2-4 .  vi 

NOMENCLATURE  FOR  SECTION  5 .  ix 

1.  INTRODUCTION  .  1 

1.1  Orographic  Effects  on  Global  Climate  ....  1 

1.2  Radiative  Transfer  in  Climatology  .  3 

2.  OROGRAPHIC  EFFECTS  .  5 

2.1  The  HAIFA  Equations .  5 

2.2  Numerical  Approximation  of  HAIFA  Equations  .  9 

2.2.1  Finite  Difference  Scheme  .  9 

2.2.2  The  Advection  Scheme .  12 

2.2.3  Update  of  Other  Terms  in  the 

Vorticity  and  Energy  Equations  ....  14 

2.2„4  Solution  of  the  Poisson  Differ¬ 
ence  Equation  by  Finite  Fourier 
Transform .  15 

2.2.5  The  FFT  Solution  of  the  Poisson 

Equation  Having  Non-Rectangular 
Boundaries .  20 

2.2.6  Description  of  Poisson  Solver 

Routines .  20 

2.3  Stability  Analysis  .  23 

2.4  Boundary  Conditions  .  25 

2.5  HAIFA  Code  Description .  28 

2.5.1  Initiating  A  Calculation  .  30 

2.5.2  Major  Subroutines  in  the  Main  HAIFA 

Calculational  Loop .  31 


iii 


3SR-79S 


TABLE  OF  CONTENTS,  contd. 

Page 

3.  MODIFICATIONS  TO  HAIFA .  33 

3.1  Compressibility .  33 

3.1.1  Derivation  of  the  Differential 

Equations .  33 

3.1.2  Method  of  Numerical  Solution  .  39 

3.2  Moisture .  41 

3.2.1  Derivation  of  Equations .  41 

3.2.2  Difference  Equations  .  47 

3.3  Variable  Zoning  in  Vertical  Direction  ....  47 

3.3.1  The  Poisson  Solver .  48 

3.3.2  Vertical  Advection  .  SI 

4.  TEST  PROBLEMS .  S3 

4.1  Single  Wave .  S3 

4.2  Two  Wave  Problem  ,  ,  ,  , . ,  67 

4.3  Uniform  Velocity  .  *  ,  .  , .  75 

4.4  Inversion  Layer  I  ......  , .  84 

4.5  Tropopause  Problem  .  .....  90 

5.  RADIATION  IN  THE  EARTH'S  ATMOSPHERE .  95 

5.1  Parameter  Specification  .  105 

5.2  Boundary  Conditions  .  113 

5.2.1  Solar  Spectrum . -113 

5.2.2  Terrestrial  Radiation  Spectrum  ....  123 


iv 


3SR-795 


TABLE  OF  CONTENTS,  contd. 

Page 

5.3  Discretization  of  the  Transport  Equation 

for  Numerical  Solution  .  125 

5.3.1  v-Discretization  .  125 

5.3.2  y-Discretization  .  135 

5.3.3  z-Discretization  .  136 

5.3.4  Numerical  Solution  .  138 

6.  FUTURE  STUDIES . 142 

6.1  HAIFA  Code  Development  and  Applications  .  .  .  142 

REFERENCES . 144 

APPENDICES 

A  -  Derivation  of  Boussinesq  Equations  .  A1 

B  -  Water  Production  Term .  B1 

C  —  Edit  Quantities  .  Cl 

D  -  HAIFA  Code  Listing .  D1 

E  —  A  Computer  Code  for  the  One-Dimensional 

Boundary  Layer  .  El 


v 


3SR-795 


a 


D 

d 

Jt 

n 


F 

g 

r 

h 

i 

J 


p 

♦ 

* 


NOMENCLATURE  for  SECTIONS  2-4 

uAt 

Ax 

specific  heat  at  constant  pressure 
drag  force  on  the  obstacle 

ft  + 

fluid  vorticity  -  -  |^- 

advective  flux  across  a  boundary 

acceleration  of  gravity 

dry  adiabatic  lapse  rate  -  g/Cp 

enthalpy 

maximum  value  of  the  grid  indice  i 
maximum  value  of  the  grid  indice  j 
numerical  grid  indices  ' 
temperature  diffusion  constant 

viscous  diffusion  constant 

latent  heat  of  vaporization  for  water 
cloud  water  content 

rain  water  content 

water  production  terms 

pressure 

compressibility  stream  function  defined  in 
Eq.  (3.14) 

stream  function 
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Q  ■  total  water  content 

q  *  water  contained  as  cloud  moisture  and  vapor 
R  **  gas  constant  for  air 
r  *  relative  humidity 

p  ■  density 

S  «  static  stahility 

s  *  entropy 

T  ■  temperature 

t  ■  time 

Vt  ■  terminal  velocity  of  water  droplet  in  atmosphere 

iu  +  kw  ■  total  velocity 

u  -  total  horizontal  velocity 
w  ■  vertical  velocity 

x  ■  horizontal  Cartesian  coordinate 
z  ■  vertical  Cartesian  coordinate 

C  ■  compressibility  vorticity  function 

■  5TC',u)  '  ^fCpw) 
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SUPERSCRIPTS 

'  =  perturbation  quantity 
"  =  defined  by  Eq.  (3.24) 
n  =  time  step  index 


SUBSCRIPTS 

D  *  diameter  of  water  droplet 

i,j  *=  numerical  grid  indices 

o  *  initial  spatial  distribution  (sometimes  used 
to  indicate  a  ground  level  value) 
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NOMENCLATURE  for  SECTION  5 

p  ■  density 
P  »  pressure 
T  ■  temperature 
z  =  vertical  coordinate 
v  ■  frequency 
X  ■  wavelength 
■  unit  vector 

6,<f>  -  spherical  angular  coordinates  defining  ft 

p  *  cos0 

Iv  »  specific  intensity  of  radiation 

Jv  »  radiation  source  function 

kv  -  volume  extinction  coefficient 

ctv  -  volume  absorption  coefficient 

o'  ■  a.  corrected  for  stimulated  emission 
v  v 

0V  ■  volume  scattering  coefficient 

Bv  ■  Planck  function 

Pv  =  scattering  phase  function 

h  *  Planck's  constant 
k  -  Boltzmann's  constant 
0g  ■  scattering  angle 

u  -  COS0e 
s  s 
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ct,v 
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a3,v 

CP 
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v  ,M 
Pn 


v  ,M 

ns 
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n(a) 


NOMENCLATURE  for  SECTION  5,  contd. 

azimuthal  average  of  1^ 

azimuthal  average  of  P^ 

radiation  energy  density 

a-component  of  radiation  energy  flux 

a3-component  of  radiation  pressure  tensor 

specific  heat  of  air  at  constant  pressure 

frequency-integrated  vertical  flux 
Rayleigh  volume  scattering  coefficient 

Mie  volume  scattering  coefficient 

Rayleigh  phase  function 

Mie  phase  function 

index  of  refraction  of  air  at  760  mm  Hg  and  15°C 

number  density  of  air  molecules 

number  density  of  air  molecules  at  760  mm  Hg  and 
15°C 

depolarization  factor  for  Rayleigh  scattering 

Mie  scattering  functions 

radius  of  (spherical)  aerosol  particle 
complex  ind.'x  of  refraction 
probability  distribution  of  aerosol  radii 
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NOMENCLATURE  for  SECTION  5,  contd. 


number  density  of  aerosols 

extinction  cross  section  of  a  spherical  particle 

scattering  cross  section  of  a  spherical  particle 

absorption  cross  section  of  a  spherical  particle 

solar  intensity 

directional -hemispherical  reflectivity 
bidirectional  reflectivity 
azimuthal  average  of 
directional  emissivity 


ground  (surface)  temperature 


ysolar 

v 


vTliff 

v 


solar-beam  part  of  I 


diffuse  part  of  Iv 
frequency -averaged  I 


Jj(vi  +  v^+1)  ■  center  of  frequency  interval 


Pi*Bi,Pi  "  Bv’Bv»Pv  evaluated  at  vi 


transmission  function 


scattering  source 


frequency-averaged  scattering  source 


q(°0 
qijk  n 


moments  of 
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1.  INTRODUCTION 

The  numerical  prediction  of  the  general  circulation 
of  the  atmosphere  predates  most  of  the  other  applications  of 
high-speed  computers  to  physical  problems.  The  codes  which 
exist  at  several  major  research  centers  have  reached  levels 
of  considerable  sophistication.  These  codes  are  used  to 
solve  time-dependent  equations  describing  atmospheric  motion 
in  a  three-dimensional  representation.  Parametric  descrip¬ 
tions  are  included  to  take  into  account  the  effects  of  inso¬ 
lation,  turbulent  transport,  and  .•cisture.  For  the  applica¬ 
tion  to  short  period  forecasts  (covering  a  time  interval  of 
several  days) ,  the  physical  processes  taken  into  account  in 
the  codes  are  quite  satisfactory,  relying  on  the  kinetic  and 
internal  energy  already  in  the  atmosphere  and  depending  less 
on  the  utilization  of  the  energy  source  from  insolation. 

For  predictions  over  a  longer  period  of  time  the  pro¬ 
cesses  which  transform  the  solar  energy  into  motion  of  the 
atmosphere  are  much  more  important.  The  many  phenomena  which 
affect  transfer  of  energy  and  moisture  through  the  earth- 
ocean-atmosphere  system  are  incompletely  described.  Descrip¬ 
tions  of  the  ocean-air,  air-land,  and  land-ocean  interfaces, 
and  of  the  topographic  boundary  conditions  are  necessary  for 
a  qualitatively  correct  predictive  model. 

1.1  OROGRAPHIC  EFFECTS  ON  GLOBAL  CLIMATE 

Phenomena  taking  place  on  a  scale  smaller  than  the 
resolution  of  global  circulation  codes  can  cause  changes  in 
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climate.  The  tropospheric  transport  coefficients  that  are 
required  in  the  global  atmospheric  model  may  arise  from  at¬ 
mospheric  motions  that  occur  in  quite  small  regions  (e.g., 
mountain  lee  waves).  Transport  is  also  effected  by  convec¬ 
tive  eddies  such  as  cumulus  and  cumulo-nimbus  convective 
cells.  These  may  be  influenced  by  small  geographic  features 
such  as  islands  and  by  upper  atmospheric  phenomena  such  as 
jet  streams  and  waves. 

The  simplest  method  of  accounting  for  meso-scale  phe¬ 
nomena  is  to  calculate  parameters  (such  as  eddy  diffusivities) 
according  to  some  fit  to  experimental  data,  risking  large 
inaccuracies  due  to  incomplete  and  inappropriate  data.  A 
technique  which  can  give  more  accuracy  is  to  compute  these 
parameters  by  means  of  several  meso-scale  calculations  per¬ 
formed  separately,  or  concurrently  with  the  large  scale  cal¬ 
culation.  This  permits  a  more  complete  description  of  rele¬ 
vant  physical  processes  to  be  built  into  the  global  model. 

Present  research  at  Systems,  Science  and  Software  (S*) 
concerns  the  development  of  a  meso-scale  code  capable  of 
studying  these  phenomena  and  in  presenting  calculational  re¬ 
sults  which  may  be  incorporated  into  global  scale  codes.  The 
basic  code,  as  discussed  in  this  report,  is  a  two-dimensional 
time-dependent  code  which  makes  use  of  the  Boussinesq  approx¬ 
imation.  The  code  is  described  in  Section  2.  Several  test 
calculations  have  been  completed  which  show  the  transient 
effects  on  the  air  flow  over  mountain  ranges  under  various 
atmospheric  conditions.  These  results  are  described  in 
Section  4.  The  momentum  transport  from  the  atmosphere  to  the 
earth  is  calculated  for  two  cases  and  these  results  are  also 
discussed. 

Modifications  to  the  code,  reported  in  Section  3, 
include  the  effects  of  moisture,  variable  zoning  i r.  the 
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vertical  direction  in  order  to  better  describe  the  atmospheric 
conditions  (i.e.,  inversion  layers,  etc.)  and  consideration  of 
the  compressibility  effects  of  the  atmosphere.  The  results  of 
test  calculations  using  these  new  codes  will  be  reported  in 
the  final  report  of  this  contract. 

1.2  RADIATIVE  TRANSFER  IN  CLIMATOLOGY 

To  quantify  the  sources  and  sinks  of  energy  in  the 
atmosphere  due  to  solar  and  terrestrial  radiation  as  a  func¬ 
tion  of  location,  season,  and  time  is  a  central  problem  in 
predictive  climatology.  Radiation  is  the  source  which  strong¬ 
ly  influences  the  level  of  response  for  all  other  parts  of 
the  system.  A  number  of  parameters  depend  sensitively  on  the 
solar  radiation:  humidity,  cloudiness,  extent  of  snow  and 
ice,  etc.,  and,  in  turn,  the  amount  of  solar  radiation  heat¬ 
ing  the  air  and  land  depends  on  the**. 

Because  of  the  intrinsic  difficulty  of  the  radiative 
transfer  calculation,  very  substantial  approximations  have 
been  made  in  the  descriptions  of  radiative  effects  in  all 
atmospheric  codes.  Most  calculations  of  atmospheric  radiation 
have  been  limited  to  approximations  of  long-wave  cooling. 

Only  a  few  transient  calculations  have  been  performed  and  the 
radiative  response  of  the  lower  boundary  of  the  atmosphere  has 
been  very  crudely  approximated  or  ignored. 

The  development  to  date  of  a  one-dimensional  radiative 
transfer  code  which  will  take  into  account  the  time-dependent 
modifications  in  the  thermal  stratification  of  the  atmosphere 
is  described  in  Section  5.  A  one-dimensional  boundary  layer 
code  which  describes  the  basic  hydrodynamics  of  the  air  flow 
and  which  will  be  used  to  test  the  radiative  transfer  code  is 
also  described.  The  radiative  transfer  code  will  be  capable 
of  characterizing  the  transfer  through  an  atmosphere  described 
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by  temperature,  pressure,  humidity,  C(K  concentration,  Oj  con* 
centration,  and  concentrations  of  other  trace  constituents 
including  aerosols. 

In  summary,  the  two  major  areas  under  investigation 
are  (1)  the  effects  of  mountain  ranges  on  energy  and  momentum 
transfer,  and  (2)  the  transient  interaction  of  solar  radiation 
with  the  earth's  atmosphere.  The  development  of  numerical 
models  to  study  these  phenomena  is  described  in  Sections  2, 

3,  and  5  of  this  report.  A  sixth  section,  describing  addi¬ 
tional  objectives  to  be  pursued  in  connection  with  these  in¬ 
vestigations  is  also  included. 
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2.  OROGRAPHIC  EFFECTS 

The  effects  of  mountain  ranges  on  atmospheric  trans¬ 
port  is  being  investigated  using  a  two-dimensional  numerical 
code  HAIFA  (hydrodynamics  in  an  Almost  Incompressible  Flow 
Approximation) .  This  code  calculates  time-dependent  dynamic 
flow  based  on  the  Boussinesq  approximation.  The  description 
of  the  basic  code  is  contained  in  this  section.  The  modifi¬ 
cations  to  this  code  to  incorporate  the  effects  of  moisture 
and  compressibility  are  discussed  in  Section  3. 

2.1  THE  HAIFA  EQUATIONS 

The  numerical  investigation  of  mountain  waves  requires 
that  the  effects  of  inertia  and  buoyancy  be  taken  into  account. 
The  two-dimensional  time-dependent  Boussinesq  equations,  devel¬ 
oped  herein,  include  these  effects  in  the  HAIFA  computer  code. 
The  buoyancy  effects  are  due  to  adiabatic  changes  of  tempera¬ 
ture  induced  by  perturbations  of  an  initially  thermally  stra¬ 
tified  atmosphere.  Deviations  from  constancy  of  the  density 
in  other  terms  of  the  fluid  equations,  including  the  continuity 
equation,  arc  neglected,  giving  a  set  of  equations  which  are 
basically  valid  for  an  incompressible  fluid.  The  use  of  the 
Boussinesq  equations  for  the  investigation  of  mountain  waves, 
therefore,  is  appropriate  in  that  the  effects  of  buoyant  sta¬ 
bility  is  restricted  by  the  incompressibility  of  the  flow. 

These  equations,  as  used  in  HAIFA,  are  the  vorticity  equation 
derived  from  the  two-dimensional  equations  of  motion ,  the 
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energy  equation,  and  the  continuity  equation  for  an  incom¬ 
pressible  fluid.  An  outline  of  the  derivation  of  these 
equations  follows.  (The  symbols  used  in  the  equations  are 
defined  in  the  Nomenclature  list.) 

In  the  Boussinesq  approximation,  the  momentum  equa¬ 
tions  in  the  horizontal  (x)  and  the  vertical  (z)  directions 
are : 


at-  *  I; £*  V,<VU>  •  t2-1) 

al*  *♦  wv*)  •  (?  2) 

For  the  present,  we  have  neglected  the  Coriolis  terms  in  this 
set  of  equations. 

The  incompressible  continuity  equation  in  two  dimen¬ 
sions  is 


3u 

3x 


♦ 


3w 

37 


0 


(2.3) 


The  vorticity  equation  used  in  the  HAIFA  code  is  derived  using 
Eqs.  (2.1),  (2.2),  and  (2.3).  Eq.  (2.1)  is  differentiated  with 
respect  to  z  and  Eq.  (2.3)  with  respect  to  x  .  Consistent 
with  the  Boussinesq  approximation,  the  variation  of  pQ  with 
height  is  assumed  negligible.  Subtracting  one  from  the  other 
removes  the  pressure  terms.  If  one  also  treats  the  diffusion 
coefficient  ky  as  a  constant,  the  resulting  expression  is: 

jk(n)  "  *  p“Q  3x  +  kv7*(n)  *  (2*4) 
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where  n  is  defined  as  the  vorticity  component  perpendicular 
to  the  x-z  plane.  Mathematically, 


n 


3u  3w 

37  '  77 


It  is  further  possible  to  modify  Eq.  (2.4)  consistent  with  the 
Boussinesq  approximations.  The  variables  p,  T  and  p  may 
be  written  as  functions  of  their  static  values  plus  a  perturba 
tion  contribution  as  follows: 


p(x,z,t)  -  p0(z)  ♦  p*(x,z,t)  , 
T(x,z,t)  -  T0(z)  ♦  T'(x,z,t)  , 
p(x,z,t)  -  p0(z)  ♦  P * (x ,  z ,  t )  . 


The  buoyancy  term 


can  then  be  written  as 


1_|£ 

P0  3x 


1_  |£ 1 
p_  ax 


(2.5) 


However,  for  the  Boussinesq  approximation  to  be  valid,  the 
density  variation  p'  must  depend  mainly  on  temperature, 
i.e.,  the  variation  of  density  due  to  t*u?  dynamical  pressure 
is  assumed  negligible  (see  Appendix  A) *  Therefore, 
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T  T 

•  • 


P'  ■  (|$)  T'  -  -  ^  T'  .  (2.6) 

Substituting  Eq.  (2.6)  into  Eq.  (2.4)  and  using  Eq.  (2.3)  to 
allow  the  result  to  be  written  in  conservative  form,  the 
vorticity  equation  is 

|f(n)  +  §^(un)  ♦  |y(wn)  “  "  f-  +  kv^2n  •  (2.7) 

o 

Eq.  (2.7)  is  the  first  of  three  equations  to  be  solved 
in  the  HAIFA  code.  The  second  equation  results  from  the  con* 
tinuity  equation  and  the  definition  of  vorticity.  Defining  a 
stream  function  such  that  u  -  3<>/3z  and  w  -  -3i|i/3x,  the 
continuity  equation  is  automatically  satisfied.  Further,  the 
stream  function  is  related  to  the  vorticity  through  a  Poisson 
equation  of  the  form 


-  n  .  (2.8) 

The  final  equation  necessary  to  complete  the  descrip¬ 
tion  of  mountain  waves  is  the  energy  equation.  This  equation 
expresses  the  first  law  of  thermodynamics 


for  an  adiabatic  system.  For  a  perfect  gas  with  constant 
specific  heat  and  using  the  hydrostatic  approximation  in  the 
dp/dt  term,  this  equation  may  be  expressed  by 


dT 

dt 


m 


(2.9) 


•  ► 
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Substituting  Eq.  (2.5)  into  Eq.  (2.9),  the  resulting  energy 
equation  is 


If-  *  ♦  sit"*) 


3T 

-w  +  F  ♦  V‘T’ 


(2.10) 


Eqs.  (2.7),  (2.8)  and  (2.10)  constitute  the  fluid  flow 
equations  integrated  in  the  HAIFA  code. 

2.2  NUMERICAL  APPROXIMATION  OF  HAIFA  EQUATIONS 

Eqs.  (2.7),  (2.8)  and  (2.10)  are  written  in  finite 
difference  form  and  integrated  numerically.  The  integration 
is  accomplished  by  updating  the  equations  in  time  for  each 
variable  based  on  the  values  at  the  previous  time  step  or  an 
intermediate  time  located  between  two  successive  time  steps. 
Each  of  these  steps  will  be  discussed  in  turn  in  this  report. 
These  descriptions  include  the  definition  of  the  grid  used  and 
the  location  of  each  variable  listed  in  the  equations,  the 
evaluation  of  the  advection  terms  in  the  vorticity  and  energy 
equations,  the  solution  for  the  stream  function  from  the 
Poisson  eq.  (2.8),  and  a  discussion  of  the  boundary  conditions 
used  in  the  numerical  integration. 

2.2.1  Finite  Difference  Scheme 

The  basic  scheme  used  to  numerically  integrate  the 
HAIFA  equations  is  shown  in  Figure  2.1.  The  finite  difference 
grid  used  in  HAIFA  is  shown  in  Figure  2.2.  The  locations  of 
the  major  variables  with  respect  to  the  grid  cells  are  defined 
in  the  figure. 
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Figure  7.2  —  HAIFA  Finite  Difference  Grid. 

The  stream  functions  are  located  at  the  grid  points,  the  vor- 
ticities  and  temperatures  are  cell  centered  and  velocities  are 
centered  on  a  grid  line  located  between  stream  line  values. 

In  this  way,  the  velocities  defined  in  finite  difference  form 
are : 


v 


ij 


_  3^  _ 

'  *ii 

(2.11) 

Yz 

Az 

2!  ... 
3x 

.  ^i.j  -  hi 

Ax 

(2.12) 
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2.2.2  The  Advection  Scheme 

The  advection  of  temperature  and  vorticity  in  HAIFA 
is  calculated  using  either  the  second  or  fourth  order  scheme 
of  Crowley. ^  The  selection  of  the  second  or  fourth  order 
scheme  is  optional  and  is  determined  by  the  trade-off  between 
accuracy  and  computing  time.  The  schemes  chosen  are  written 
in  conservation  form  and  are  based  on  forward  time  differences 
and  centered  space  differences.  Test  calculations  performed 
by  Crowley  indicated  that  for  the  same  order  of  accuracy,  the 
conservation  form  produced  more  accurate  solutions  than  the 
advection  form. 

In  the  conservation  form,  the  time  derivative  and  ad¬ 
vection  terms  of  the  vorticity  or  temperature  equation  may  be 
written  as 


.  |IU£I  .  |Mi  .  s 


(2.13) 


where  4)  is  either  T  '  or  n  and  S  is  the  source  term. 

In  two  dimensions  a  splitting  technique  is  used;  the 
calculational  scheme  calls  for  solving  a  one-dimensional  equa¬ 
tion  twice,  i.e.,  the  net  flux  of  vorticity  or  temperature  is 
solved  for  in  the  horizontal,  the  quantity  solved  for  in  the 
zone  being  updated  due  to  this  flux  and  the  procedure  is  then 
repeated  in  the  vertical  direction  using  the  partially  updated 
values.  The  equation  for  the  flux  across  the  boundary  j 
written  in  finite  difference  form  (second  order  accurate)  is 


) 


where 


Uj  At/Ax 


(2,14) 
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'he  net  change  in  the  variable  4  in  the  cell  ij  due  to  ad* 
vection  in  the  horizontal  is  then 

♦ntl  -  ♦"  -  j|(Fj  -  FjM)  .  (2. IS) 

I 

The  corresponding  fourth  order  scheme  for  the  flux 
across  the  boundary  j  is 

Fj  “  *  ♦j-l5  ‘  (*j+l  *  *j-2J] 

a  2 

*  *  ♦j-P  ‘  (*j+l  *  ^j-2^] 

a  I 

*  t4~  [(*J  *  ♦j-l5  *  (Vl  ■  ♦j-2)] 

+  -  ♦)-!>  •  <♦)♦!  *  ♦j-2']J  •  C2-16> 

The  numerical  stability  of  these  equations  is  discussed  in 
Section  2.3  of  this  report.  The  accuracy ,  as  discussed  by 
Crowley,  is  found  by  expanding  the  quantities  in  Taylor  series, 
both  in  time  and  in  space.  The  result  gives  the  solution  of 
the  variable  $  at  the  new  time  accurate  to  order  At2 
in  time.  The  time  derivative  of  the  finite  difference  form 
of  the  differential  equation  is  thus  accurate  to  order  At* 
in  time.  The  second  order  scheme,  Eq.  (2.14),  has  a  trun¬ 
cation  of  order  Ax3  and  the  fourth  order  scheme,  Eq.  (2.16), 
is  accurate  to  Axs  in  space. 
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2.2.3  Update  of  Other  Terms  in  the  Vorticity  and  Energy 
Equations 

The  vorticity  equation  has  two  additional  terms  besides 
the  advective  terms.  In  general,  central  differences  are  used 
in  the  numerical  scheme.  The  buoyancy  term 


is  expressed  as 


*  T?“ 


'  lEx 


The  diffusion  term  k  V*n  is  expressed  as 


kv|(ni,j*l  "  2nij  * 


r  At 


*  (ni*l,j  "  2nij  ^  ni-l,j)(2x*) 


(2.17) 


(2.18) 


The  energy  equation  is  handled  in  a  similar  manner. 
The  diffusion  term  ktV*T  is  expressed  as  in  Eq.  (2.18) 
with  q  replaced  by  T  .  The  remaining  term  in  the  energy 
equation 


-w 


r 


is  calculated  using  centered  quantities.  The  term  in 
brackets  is  calculated  analytically  from  input  data  and  the 
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I 

I 

r 

r 

i 

r 

i 

i 


velocity  is  expressed  as  an  averaged  quantity 


w 


*  Wij 


2.2.4  Solution  of  the  Poisson  Difference  Equation  by  Finite 

Fourier  Transform 

The  solution  of  the  Poisson  equation  by  neans  of 
Fourier  transform  results  in  a  direct  (or  exact)  solution  of 
the  difference  equations  and  their  boundary  values.  In  the 
current  version  of  the  subroutine  there  are  some  limitations 
on  the  generality  of  the  solution;  the  spatial  interval  Ax 
must  be  constant  (see  Section  3.3  for  variable  Az).  The  solu¬ 
tion  must  be  periodic  in  the  x-direction  and  prescribed  values 
of  the  stream  function  are  to  be  maintained  on  the  top  and 
bottom  boundaries  of  the  rectangular  region.  Hon  the  hounds 
ary  is  modified  from  the  rectangular  shape  is  discussed  in 
the  following  section.  • 

A  second  order  finite  difference  approximation  to  the 
Poisson  equation  V2ip  »  n  is  obtained  by  replacing  the  second 
derivative  operator  by  a  centered  second  difference  operator. 


tk  hi  •  f 1  *ii 

(Ax)1  (At) J 


i  •  1,2,  ...,  1 
i  *  2,  ...,  J-l  , 


where 


and 


-  2*i)  *  *i-i, i 

*i  *ij  *  *!,)♦!  ’  2*ij  4  *1  J-l  • 


(2.19) 
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Boundary  conditions  are  imposed  as  follows: 
At  the  bottom  of  the  mesh, 

j  .  i  ■  1 ,  ...»  I  . 

At  the  top  of  the  mesh. 


♦i  j  ■  ,  i  •  1,2,  . . . ,  I  . 

The  cyclic  boundary  conditions  in  the  horizontal  are, 

*o,j  ’  *I,j  and  *1,J  ■  *1.1, )  •  i  ’  2 . J'1  • 

We  introduce  an  orthonormal  base  set  of  functions  having 
cyclic  properties  on  the  index  i: 

wik  "  cos  ,  I  is  even 

wi  I-k  "  sin  2fkl  »  i  ■  1,2,  ...,  I  . 

wi,I  "  • 

k  ■  1,2,  (I,,  y  *  1  • 

wi,I/2  "  1 W  cos  1  » 

These  are  the  finite  Fourier  functions  which  have  the  properties, 

I 

£  wik  “a  ■  4ki 

i-1 
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and  the  analogous  cyclic  boundary  conditions  are  valid  in  the 
horizontal.  They  also  have  the  property  that  they  are 
eigenfunctions  of  the  central  second  difference  operator 


6 


i 

x 


wik 


wik 


where  A^  ■  2  sin  irk/ 1  .  These  functions  are  complete  func¬ 

tions  on  the  interval  i  -  1,2,  ...»  I  .  Consequently, 
an  arbitrary  function  f^  on  this  space  can  be  represented 

I 

£i  -  L  *k  wik 
k-l 


where 


A 

t 

i-1 


fi  Kik 


We  are  now  ready  to  consider  Eq.  (2.19)  from  the  point 
of  view  of  Fourier  transformation.  The  vorticity  and  stream 
function  are  represented  as  Fourier  series  as  follows; 


and 


& 

£ 


i 

r 

k-l 

I 


bkj  wik 


2  ak;  wik 
k-l 


where 


ik 


where 


I 


2  *ij 

i-  1 


wik 


(2.20) 
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Substituting  into  Eq. 


(2.19)  we  obtain 


Multiplying  by  w^  and  summing  over  i  gives 


j  "  2 ,  •  •  •  ,  *1*1  p 
i  -  1,2,  .  .  * ,  I  . 


(2.21) 


The  values  of  a1>t  and  ajt  required  by  Eq.  (2.21)  are 
obtained  from  the  boundary  values 


I 

•i.»  ■  £  wit  and 
i-1 

l  (2.22) 

•j,t  ■  X  bi  wu  • 

i-1 


In  Eq.  (2.21)  the  value  of  the  wave  number,  1,  appears 
only  parametrically.  For  each  value  of  l  there  is  a  tri¬ 
diagonal  equation  having  fixed  values  at  the  end  points  of  the 
j -interval. 

We  summarize  the  procedure  for  obtaining  the  direct 
solution  of  the  Poisson  equation,  Eq.  (2.19),  by  Fourier  trans 
form: 


18 


3SR-795 


(1)  The  vorticity  and  the  top  and  bottom  boundary 
values  of  the  stream  function  are  subjected  to  Fourier  trans¬ 
formation  to  obtain 


I 

bjt  "  nij  wik  » 
i-1 

1 

•l,i  *  £  °i  Wii  •  and 
i-1 

I 

aJ,i  "  2  6i  WU 
i-1 

(2)  The  Fourier  components  of  the  stream  function  are 
obtained  by  solving  the  tridiagonal  system  of  equations,  Eq. 
(2.21),  for  a^. 

(3)  The  stream  function  itself  is  obtained  by  Fourier 
synthesis 

I 

^i j  "  aj  £  wi£ 

l-l 

The  quantity  I  must  be  even.  In  order  to  take  maxi¬ 
mum  advantage  of  the  efficiency  of  the  Fast  Fourier  Transform, 
the  quantity  I  should  also  be  a  power  of  2. 
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2.2.5  The  FFT  Solution  of  the  Poisson  Equation  Having  Non- 

Rectangular  Boundaries 

In  order  to  represent  a  mountain  within  the  computa¬ 
tional  grid  it  is  necessary  to  depart  from  rectangular  bound¬ 
aries.  A  modification  of  the  solution  algorithm  using  the  FFT 
is  necessary  to  take  account  of  the  specified  values  of  i|;  on 
the  mountain  contour.  The  procedure  for  carrying  out  this 

modification  of  the  direct  solution  of  Poisson's  equation  on 

(Z) 

an  irregular  region  has  been  described  by  Buzbee,  et.al.v  1 

We  consider  the  case  in  which  there  are  p  internal 
grid  points  on  which  the  potential  is  to  be  specified.  These 
points  constitute  the  adjacent  mesh  points  lying  along  the 
boundary  of  the  mountain  which  will  be  assigned  the  same  value 
of  potential  (usually  zero)  as  the  lower  boundary.  The  first 
step  is  to  precalculate  the  stream  function  contribution  at 
each  of  the  p  points  of  unit  vorticity  located  at  each  of 
the  points.  The  solution  is  then  obtained  ,  solving  Poisson's 
equation  twice  for  each  cycle.  First,  Poisson's  equation  is 
solved  with  arbitrary  vorticity  on  the  boundary  points.  The 
difference  between  the  obtained  and  desired  values  of  the 
stream  function  at  e.  zh  of  the  p  points  is  used  to  obtain 
the  corresponding  vorticity  increments  through  application  of 
the  precalculated  matrix.  A  second  solution  of  Poisson's 
equation  using  the  incremented  vorticity  field  gives  the  final 
value  of  the  stream  function  within  the  calculational  region. 

2.2.6  Description  of  Poisson  Solver  Routines 

This  section  describes  the  subroutines  currently  used 
in  the  HAIFA  code  to  solve  the  Poisson  equation  in  x-z  geometry. 
The  method  of  solution  employs  a  Fourier  transform  in  the 
x-direction,  solving  the  resultant  set  of  one-dimensional  dif¬ 
ference  equations  (one  for  each  wave  number)  by  Gaussian 
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elimination  in  the  z-direction  and  performing  the  inverse 
x-direction  Fourier  transform  to  obtain  the  solution.  The 

(3) 

Cooley -Tukey  Fast  Fourier  Transform  (FFT)  technique 
is  employed  (subroutine  COOTUK)  with  some  pre-  and  post¬ 
processing  of  the  data  for  efficient  utilization  of  the  al¬ 
gorithm.  In  the  current  version  the  dependent  variable  (the 
stream  function  \p  in  the  HAIFA  context)  is  assumed  to  have 
cyclic  boundary  conditions  in  the  x-direction  and  fixed  values 
at  che  top  and  bottom  of  the  grid. 

At  the  beginning  of  each  new  calculation,  there  are 
references  to  subroutines  which  are  used  only  once  in  each 
problem.  These  are  called  SETUP  and  OBSET. 

SETUP  --  This  entry  references  an  internal  subroutine  SET, 
whose  function  is  to  define  certain  index  parameters  and  re¬ 
quired  data  arrays  that  are  used  throughout  the  calculation 
by  the  Poisson  solver, 

OBSET  --  This  subroutine  is  called  only  when  internal  boundary 
conditions  are  to  be  applied.  Suppose  there  are  p  internal 
points  required  to  have  stream  function  values  ij>° ,  ,  ...  i|>p. 

This  subroutine  computes  a  p  x  p  matrix  C  which  has  the  fol¬ 
lowing  property: 

a  unit  vorticity  is  placed  in  the  posi¬ 
tion  of  internal  boundary  point  j .  The  value 
of  the  independent  variable  (vorticity)  is 
assumed  to  be  zero  at  every  other  point.  The 
Poisson  equation  solver  XYPQIS  (see  discus¬ 
sion  below)  is  called  and  returns  the  influ¬ 
ence  of  that  particular  unit  vorticity  on  all 
the  other  internal  boundary  points.  These 
influences  are  put  into  row  j  of  matrix  C. 
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it  m 


} 

L_; 


*  1 


This  procedure  is  continued  until  all  p  in¬ 
ternal  boundary  influences  have  been  computed. 
Finally,  subroutine  OBSET  forms  and  stores 
the  inverse  matrix  C~ 1 . 

The  controlling  subroutine  for  the  Poisson  equation 
solution  is  named  LAPLAC  (for  the  Laplacian  symbol  V2). 

This  routine  is  responsible  for  the  solution  to  both  standard 
boundary  condition  cases  and  problems  which  include  internal 
boundaries . 

Each  cycle,  subroutine  LAPLAC  averages  the  cell- 
centered  HAIFA  vorticities  to  provide  node-centered  vortici- 
ties.  Then  the  Poisson  equation  solver  XYPOIS  is  called  to 
provide  the  updated  values  of  the  stream  function.  In  the 
case  of  internal  boundaries,  one  more  step  is  performed  in 
subroutine  LAPLAC.  Upon  the  first  return  from  solving  the 
Poisson  equation,  each  internal  boundary  has  a  value 
i-1,  ...  ,p  which  in  general  is  not  the  required  value  i|»?. 

A  ector  A^J>  of  the  differences  is  formed.  Then, 

using  the  inverse  matrix  C’1  formed  in  subroutine  OBSET,  one 
may  compute  the  required  modifications  Aq^  to  the  values  of 
the  independent  variable  at  each  of  the  p  internal  boundary 
points  from 


The  independent  variable  is  so  modified,  and  the  XYPOIS 
package  is  called  once  again.  The  solution  returned  now  con¬ 
tains  the  correct  values  for  the  internal  boundary  points  as 
well  as  the  other  grid  points.  It  remains  to  discuss  the  sub¬ 
routine  XYPOIS. 
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XYPOIS  --  This  entry  is  used  every  calculational  cycle  to 
carry  out  the  solution  of  Poisson's  equation.  It  contains 
as  an  argument  the  values  of  the  inhomogeneous  term  (here, 
vorticity)  in  the  interior  (nodal)  points  of  the  grid,  and 
the  fixed  values  of  the  dependent  variable  (here,  the 
stream  function)  at  the  top  and  bottom  of  the  grid.  XYPOIS 
references  four  internal  subroutines: 

(1)  FFANL  (fast  Fourier  analyzer),  which  is  respon¬ 
sible  for  carrying  out  the  x-direction  transform  of  vorticity 
into  Fourier  components.  It  processes  two  rows  at  a  time, 

so  an  uncoupling  of  the  row  components  is  required  upon  re¬ 
turn  from  the  FFT  routine  COOTUK; 

(2)  GAUSS,  which  is  responsible  for  solving  the  re¬ 
sulting  z-direction  tridiagonal  equations  for  the  transform 
of  the  dependent  variable  (see  Section  3.3.1); 

(3)  FFSYN  (fast  Fourier  synthesizer),  which  is  the 
inverse  of  FFANL,  is  responsible  for  restoring  the  Fourier 
components  to  the  new  values  of  the  independent  variable  by 
another  call  to  subroutine  COOTUK.  These  values,  represents 
ing  the  solution  to  the  Poisson  equation,  are  returned  to 
the  calling  routine  (subroutine  LAPLAC)  in  the  array  contain¬ 
ing  the  original  argument  list;  and 

(4)  COOTUK,  which  carries  out  the  Cooley-Tukey  fast 
Fourier  transform. 

2,3  STABILITY  ANALYSIS 

A  numerical  stability  analysis  of  the  advection  terms 

in  the  vorticity  and  temperature  equations  has  been  completed 

r i  4} 

by  other  researchers.  Among  them,  Crowley'  *  J  did  a  complete 
analysis  for  the  scheme  presently  being  used  in  the  HAIFA 
code.  The  results  obtained  by  Crowley  indicate  that  both  his 
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second  and  fourth  order  scheme  are  stable  for  all  wave  num 
bers  if 


uAt 

Ix~ 


<  1  . 


Further,  the  fourth  order  conservation  scheme  being  used  in 
HAIFA  is  stable  for  (uAt/Ax)  <  l.S. 

As  indicated  by  Crowley,  the  schemes  both  result  in 
amplitude  drvping  and  phase  lag.  For  long  wavelength  dis¬ 
turbances  the  damping  and  phase  errors  are  appreciably  smaller 
for  the  fourth  order  scheme  than  for  the  second  order.  Com¬ 
parison  tests  with  a  typical  mountain  wave  problem  indicated, 
however,  that  the  differences  between  fourth  and  second  order 
solutions  are  not  large.  Most  of  our  calculations  have  been 
performed  with  the  second  order  scheme.  The  criterion 
built  into  the  HAIFA  code  is  more  stringent  than  any  of  those 
noted  above ,  i .e . , 


uAt 

Ix~ 


<  0.8 


A  stability  criterion  also  has  been  established  for 
the  diffusion  terms,  however,  in  all  problems  calculated  for 
this  research,  the  diffusion  coefficients  are  set  to  zero  and 
thus  these  terms  play  no  part  in  the  solution. 

One  unstable  region  was  found  using  the  above  cri¬ 
teria  in  computing  the  uniform  velocity  problem  discussed  in 
Section  4.3.  The  details  of  the  instability  and  the  new 
criteria  developed  for  that  problem  are  also  given  in  that 
section. 
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2.4  BOUNDARY  CONDITIONS 

The  initial  value  problem  solved  using  the  HAIFA  code 
requires  initial  temperature,  vorticity  and  stream  function 
distributions.  This  is  accomplished  by  prescribing  a  value 
of  the  stream  function  which  is  constant  in  the  horizontal 
direction  and  which  gives  the  desired  horizontal  velocity 
distribution  as  a  function  of  the  vertical  coordinate.  The 
vertical  velocity  component  is  set  to  zero.  The  vorticity 
at  each  point  in  the  grid  is  calculated  analytically  using 
the  definition 

3u  .  3v  .  , 

n  ■  »  since  1S  everywhere  zero. 

The  temperature  distribution  is  specified  as  being  horizontal - 
ly  stratified  with  a  lapse  rate  which  may  vary  with  altitude. 
It  is  also  possible  to  simulate  inversions. 

At  the  beginning  of  the  calculation,  with  the  flow 
already  established,  an  obstacle  is  placed  in  the  stream  by 
setting  the  lower  surface  streamline  to  coincide  with  the 
mountain  surface.  A  rigid  lid  (constant  streamline)  is  im¬ 
posed  on  the  upper  boundary  of  the  problem.  Figure  2.3  indi¬ 
cates  these  boundary  conditions  in  graphical  form. 

The  boundary  condition  imposed  at  the  sides  of  the 
grid  assumes  the  flow  to  be  cyclic,  i.e.,  the  stream  function 
at  each  vertical  grid  line  j  on  the  left  side  of  the  grid 
is  set  equal  to  the  corresponding  stream  function  at  the 
right  side  of  the  grid.  Mathematically ,  this  can  be  ex¬ 
pressed  as  41^  *  j  •  A  graphical  explanation  of  this 

boundary  condition  is  also  given  in  Figure  2.3. 
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♦  •  constant 


Upper  and  Lower  _  _ 
Boundary 

Condition  - 


■  constant 
(usually 
zero) 


Cyclic  Boundary 
Condition 


Figure  2.3  —  Schematic  of  HAIFA  Boundary  Conditions. 
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One  further  boundary  condition  is  necessary  to  obtain 
the  transient  solution.  The  vorticity  equation  requires  that 
the  temperature  gradient  in  the  x-direction  be  specified  at 
the  cell  center  bounded  by  the  obstacle.  This  requires  a 
value  for  the  temperature  perturbation  on  the  obstacle 
boundary.  The  assumption  is  made  that  the  air  immediately 
next  to  the  mountain  has  risen  from  the  bottom  of  the  grid. 
The  temperature  of  the  air  alongside  the  mountain  is  thus 
given  by 

T  -  T  -  r  •  z. 

m  o 

where 


T  -  f.he  temperature  along  the  vertical 
mountain  boundaries 

Tq  *  the  temperature  at  ground  level 

r  «  the  dry  adiabatic  lapse  rate 
z  -  the  distance  above  ground  level. 


Since  the  initial  temperature  profile  (T^)  is  given  as 
an  analytic  function  of  z,  the  temperature  perturbation  along 
the  mountain  is 

T*  -  T  -  T*z  -  T.  . 
mo  1 


Referring  to  the  example  oi  a  two-cell  thick  mountain  in  the 
figure  below,  3T'/3x  at  the  cell  centers  adjacent  to  the 
obstacle  are  calculated  as 


Tm 

n 


Ti 


T* 

4-1 


Ax 


i 
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3Ti*3 
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T '  ♦  T' 

Ti+3  *i*4 


-  T 


m 


"SF 


i+1  i+2  i+3  i+4 


At  cycle  zero  (time  equal  to  zero) ,  these  boundary  con 
ditions  are  then  used  to  determine  the  new  distribution  of 
streamlines  within  the  calculational  grid.  This  completes  the 
required  information  to  start  the  computation. 

2.S  HAIFA  CODE  DESCRIPTION 

A  flow  chart  giving  the  calculational  sequence  of  the 
HAIFA  code  is  displayed  in  Figure  2,4.  A  description  of  how 
problems  are  generated  and  the  major  subroutines  within  the 
code  is  presented  below. 
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READ  INITIAL 
INPUT  DATA 


(INPUT) 


(RTAPE) 


(OBSET) 


(EDIT) 

(OUTPUT) 

(WTAPE) 


Figure  2.4 


Flow  Diagram  of  HAIFA  Code. 
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2.S.1  Initiating  A  Calculation 

There  are  two  methods  for  initiating  a  calculation: 
generating  a  new  problem,  and  restarting  a  partially  completed 
calculation  from  a  data  tape.  These  are  controlled  by  sub¬ 
routine  INPUT. 

Generating  a  New  Problem  -  Suoroutine  INPUT  reads  all 
input  data  and  sets  up  several  constants  which  will  be  used 
in  the  calculation.  The  initial  streamline  distribution  is 
computed  from  a  series  of  input  parameters,  MT  ,  DTI  ,  DT2  , 
DT3  ,  DT4  ,  and  ZETA  such  that 

*(z)  -  DTI  ♦  DT2«zMT  ♦  DTS^^*1) 

♦  DT4*exp(-ZETA*z)  . 

These  parameters  define  the  horizontal  velocity  distribu¬ 
tion 


u(z)  -  -  MT-DT2-z(MT'1)  +  (MT+l)*DT3*zMT 

-  DT4 • ZETA* exp (- ZETA *z)  . 

The  initial  vorticities  are  found  from  differentiating  the 
above  expression  with  respect  to  z  ,  i.e.,  n  *  3u/3z  since 
3v/3x  is  everywhere  zero  at  time  equal  zero. 

The  initial  temperature  distribution  is  set  in  a 
similar  fashion  using  the  input  parameters  KT  ,  ATI  ,  AT2  , 
AT 3  ,  AT 4  ,  and  ALPHA. 

T(z)  «  ATI  +  AT2*zKT  +  AT3«z(KT+1) 

+  AT4 • exp (-ALPHA* z)  . 
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Internal  Boundaries  —  The  input  variable  NOBS  defines 
the  number  of  internal  points  which  are  to  have  fixed  stream- 
function  values.  A  series  of  data  cards  specifying  the  grid 
points  and  the  associated  values  are  read  if  NOBS  >  0  . 

Such  internal  boundary  points  are  used  to  define  grid 
obstacles,  which  are  outlined  by  a  series  of  connected  points. 
Typically,  the  fixed  value  of  ijj  assigned  to  the  obstacle 
points  is  the  lower  boundary  streamfunction  value.  The  re¬ 
quested  initialization  of  the  streamfunction,  vorticities,  and 
velocities  in  the  case  of  internal  boundaries  is  handled  by 
subroutine  OBSET. 

Restarting  A  Calculation  —  The  option  to  restart  a 
calculation  is  keyed  by  the  input  parameter  RESTRT.  If  it 
is  non-zero  in  value,  the  data  tape  is  scanned  in  subroutine 
RTAPE  until  the  cycle  requested  by  input  parameter  ISTART  is 
found.  The  values  of  the  necessary  calculational  variables 
of  the  requested  cycle  are  then  read,  and  the  computation  is 
continued. 

2.5.2  Major  Subroutines  in  the  Main  HAIFA  Calculational  Loop 

UPDATE  --  UPDATE  is  used  to  solve  the  conservative  equations 
for  vorticity  and  temperature.  Crowley's  second  order  or 
fourth  order  scheme  is  called  from  this  subroutine  to  calculate 
the  advection  terms.  This  scheme  is  described  in  Section  2.2.2 
of  this  report. 

LAP LAC  --  The  Pcisson  equation  relating  the  stream  function 
and  the  vorticity  is  solved  using  this  subroutine  as  the  con¬ 
trolling  program.  The  details  of  the  Poisson  solver  are  given 
in  Section  2.2.4. 
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VELOC  --  The  updated  stream  function  values  are  differenced 
~r  in  z-space  to  provide  the  horizontal  velocity  field  u  ,  and 

**»  in  x-space  to  provide  the  vertical  velocity  field  v. 

m.  PRTTST  --  This  subroutine  defines  the  type  of  output  required 

in  each  cycle,  viz,  plots,  large  edits,  and/or  data  dumps  on 
tape  are  available  options  with  this  program. 

TIMSTP  --  The  TIMSTP  subroutine  calculates  a  time  step  to 
be  used  in  the  calculation  limited  by  the  numerical  stability 
criterion.  The  stability  criterion  is  outlined  in  Section  2.3. 


i'j 

<*  * 
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3.  MODIFICATIONS  TO  HAIFA 

The  HAIFA  equations  described  in  the  preceeding  sec¬ 
tions  are  limited  in  that  the  formulation  has  been  simplified 
both  from  the  mathematical  and  physical  points  of  view.  In 
Section  3  we  discuss  several  investigations  to  generalize  both 
the  mathematical  and  physical  aspects  of  the  code.  The  three 
programs  described  below  are  currently  being  tested  and  are 
approaching  operational  status.  Additional  features  are  to  be 
incorporated  in  the  latter  part  of  the  contract;  they  are  dis¬ 
cussed  in  Section  6. 

3.1  COMPRESSIBILITY 

3.1.1  Derivation  of  the  Differential  Equations 

The  use  of  HAIFA  for  the  investigation  of  mountain 
waves  is  appropriate  in  that  the  effects  of  buoyant  stability 
and  dynamics  are  taken  into  account,  but  its  applicability  is 
restricted  by  the  incompressibility  of  the  flow.  In  particu¬ 
lar,  if  the  height;  of  the  mountain  range  is  comparable  with 
the  atmospheric  scale  height  there  will  be  effects  induced  by 
the  expansion  experienced  by  an  air  packet  in  being  lifted  over 
the  mountain. 

The  effects  of  compressibility  are  to  be  determined 
through  the  use  of  a  new  code  developed  with  which  problems  in¬ 
cluding  this  effect  may  be  run  and  the  results  compared  with 
HAIFA  calculations.  Several  objectives  were  sought  in  arriving 
at  a  method  of  accomplishing  this  task.  They  are  discussed 
below. 
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(1)  Sound  waves  should  be  excluded  from  the  numerical 
solutions  in  order  to  permit  efficient  calculations  having  time 
intervals  comparable  with  material  displacement  through  a  space 
interval . 

(2)  Compressibility  effects  should  be  retained. 

(3j  The  scheme  should  be  formulated  in  physical  vari¬ 
ables  to  facilitate  addition  of  new  physical  effects  (such  as 
Coriolis  force  or  water  vapor)  . 

(4)  Conservative  difference  equations  should  be  sought. 

(5)  The  scheme  should  retain  a  mathematical  form  simi¬ 
lar  to  HAIFA  to  make  programming  and  check-out  as  speedy  as 
possible . 

The  "anelastic"  equations  of  Ogura^  meet  some  of  the  above 
criteria  and  will  be  compared  further  below.  However,  the 
anelastic  equations  do  not  allow  an  arbitrary  atmospheric 
stratification,  do  not  include  the  change  in  density  due  to 
temperature  perturbations  and  are  formulated  in  problem- 
dependent  variables.  These  limitations  can  be  avoided,  as 
indicated  below. 

The  equations  for  inviscid  fluid  flow  on  a  non¬ 
rotating  earth  (additional  terms  will  be  discussed  later)  are 
written  in  conservative  form  as  follows: 


it *r-§f ♦  . 

(3.1) 

If*  *  IF1  *  IP  -  IF  *  If  - »  . 

(3.2) 

|£V  ,  |£UV  ,  |£V1  ,  |fWV  ,  |£  .  „  f 

(3.3) 

||W  +  |£UW  ,  |p  +  |fwl  ,  H  .  ,gp  _ 

(3.4) 
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For  a  perfect  gas,  ^ 
written 


RT 

P 


and  the  energy  equation  can  be 


dT  .  1  dp 

eft  "  pCp  at 


(3.5) 


For  the  adiabatic,  inviscid  non-rotating  motion  of  a  perfect 
gas  we  have  the  equations  of  motion  given  by  Eqs .  (3.1) 
through  (3.5). 

We  now  consider  the  equations  in  two  spatial  dimen¬ 
sions,  the  vertical  and  down -wind  directions,  assuming  that 
there  is  no  dependence  of  any  quantity  on  the  y- direction. 
The  equations  become 


8p  +  3£U  +  3pw  =  0 

3 1  3x  3  z  u  » 

(3.6) 

3£U  3£u1  +  3£WU  3£  .  o 

3t  3x  3  z  3x  u  * 

(3.7) 

IfiW  ,  |£UW  ,  |fwl  t  Sf  -  -gp  , 

(3.8) 

dt  ■  pc  Ht  ’  P  -  pRT  . 

(3.9) 

We  consider  the  problem  in  which  the  mountain  perturbs 
an  initially  steady  state  of  the  atmosphere,  wq  =  0, 
uQ  =  u(z  ,t=0)  ,  Tq  =  T(z , t~0)  ,  pQ  =  p(z ,t=0) ,  the  normal 
initial  conditions  previously  discussed.  These  initial  values 
are  related  to  each  other  through  the  static  atmosphere  equa¬ 
tion, 
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3P0 

ST"  *  '*po 


«>o 

ITT*  * 

o 


3  lii  p_  _ 

Ti - “  rib  • 

o 


(3.10) 


The  transient  solution  is  obtained  by  solving  the 
equations  of  motion  starting  with  the  initial  values  and  im¬ 
posing  boundary  values  on  the  motion.  In  order  to  eliminate 
sound  waves  from  the  transient  solution  it  is  sufficient  to 
set  *  0  in  Eq.  (3.6).  The  resulting  equation. 


3pu  3pw 
3x  3z 


0 


(3.11) 


is  in  suitable  divergence-free  form  for  the  introduction  of 
a  solinoidal  function.  Since  Eq .  (3.11)  is  no  longer  suffi¬ 
cient  to  determine  how  the  density  changes,  it  is  necessary 
to  introduce  an  additional  equation  based  on  an  approximation. 
We  assume  that  the  density  can  be  determined  at  every  position 
from  the  perfect  gas  equation  of  state  in  which  the  pressure 
takes  the  value  associated  with  the  static  atmosphere,  pQ  , 
through  the  relation 


P 


Po 

Rl 


(3.12) 


The  temperature  equation  can  be  written  in  terms  of 

the  deviation  T'  of  the  temperature  from  its  static  value 

(T  =  T  +  T')  . 
o' 

In  addition,  the  temperature  equation  can  be  reformu¬ 
lated  so  that  the  advection  term  assumes  a  conservative  form. 
Expanding  the  left-hand  term  of  Eq.  (3.9) 
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dT 

at 


3T ' 

TT~ 


3T ' 
u5x~ 


w3  T* 

*TT 


+  W' 


dTo 

cfz- 


f 


multiplying  by  p  and  adding  Eq.  (3.6)  multiplied  by  T'  , 
we  obtain 


3pT'  .  3pT’u  3pT'w  irdTo 

tt~  +  3^ —  +  7i —  +  war- 


1  d£ 

C  3t  * 


Assuming  the  pressure  to  have  its  static  value  in  the  right- 
hand  term  of  Eq.  (3.9), 


dp  dp0 


dP, 


ar-  *  war" "  ‘w«p< 


the  energy  equation  becomes 


3pT*  .  3pT'u  .  3pT'w 
7T~  3x - 71 - 


3T 


-w  Pg_  +  por 


(3.13) 


where  we  have  assumed 

po  T  To  +  T’ 

T  =  =  ~ r — 

o  o 

from  Eq.  (3.12) . 

The  "anelastic"  equations  of  Ogura  and  Phillips  also 
take  account  of  compressibility  effects  in  the  atmosphere  and 
it  is  of  interest  to  compare  the  above  development  with  them. 
The  anelastic  equations  are  basad  on  several  assumptions: 
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(1)  The  potential  temperature  is  almost  con¬ 
stant;  deviations  from  constancy  are 
small . 

(2)  The  density  appearing  in  Eq.  (3.11)  is 
that  associated  with  a  neutrally  strati¬ 
fied  atmosphere. 

(3)  The  potential  temperature  appearing  in 
the  momentum  equations  is  that  of  the 
neutral  atmosphere,  i.e.,  it  is  constant. 

This  assumption  corresponds  to  using  a 
neutral  atmosphere  density  in  the  ad- 
vection  terms  of  Eqs.  (3.7)  and  (3.8). 

The  treatments  of  the  buoyancy  term  of  the  momentum  equations 
and  the  energy  equations  are  the  same  in  the  two  schemes. 
Consequently,  the  proposed  scheme  is  more  general  in  two 
principal  respects;  the  initial  stratification  of  the  atmos¬ 
phere  can  be  arbitrarily  specified,  and  the  effect  of  temp¬ 
erature  changes  in  the  atmosphere  are  reflected  in  all  of  the 
density  terms. 

The  system  of  compressibility  equations,  Eqs.  (3.7), 
(3.8),  (3.11),  and  (3.13),  have  a  form  similar  to  the 
Boussinesq  equations,  and  can  be  solved  in  a  similar  manner. 
From  Eq.  (3.11)  a  stream-function-like  quantity  4>  can  be 
introduced: 


In  terms  of  a  vorticity- like  function  x,  , 

5pu  5pw 
;  '  3z  "  Tx~  * 


(3.14) 


(3.15) 
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yielding  the  same  Poisson  equation  as  for  the  Boussinesq 
approximation , 


C 


3  4» 


S3F  T  Tz 


& 


(3.16) 


The  prognostic  equation  for  C  is  obtained  by  cross  differ¬ 
entiating  Eqs .  (3.7)  and  (3.8)  and  subtracting. 


(3.17) 


Eq.  (3.17)  replaces  the  vorticity  equation  of  the 
Boussinesq  equations,  differing  principally  in  having  the 
additional  terms  containing  the  derivations  of  <}>,  u  and  w. 


3.1.2  Method  of  Numerical  Solution 

The  system  of  compressibility  equations  is  seen  to  be 
very  similar  to  the  HAIFA  equations,  and  in  fact  only  nominal 
modifications  to  the  HAIFA  code  were  required  to  produce  a 
compressible  low-speed  flow  code. 

Generating  Initial  Conditions  —  As  with  HAIFA,  the 
values  of  uQ(z)  and  Tq(z)  are  specified  by  input  to  the 
code.  In  addition,  the  initial  surface  pressure  PQ(z=0) 
must  be  specified.  The  remaining  initial  pressures  are  found 
using  the  relation  of  Eq.  (3.10), 


p0(z)  =  po(z=0)  exp  -  | 


/ 


dz 

nrr 
0  v  J 
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The  initial  density  profile  then  follows  from 

p0CO 

po(z)  =  RT  (z)  • 

o 

The  stream-function-like  quantity  <{>  is  formed  by  integrating 

■  P0CO  uo(z)  , 

and  the  vorticity-like  quantity  £  is  initialized  from  the 
non-zero  term  of  Eq.  (3.15), 


(as  in  HAIFA,  there  is  no  x-dependence  of  any  quantity  initially) . 

Modification  of  Advection  Scheme  —  The  quantities  to 
be  advected  in  the  system  of  compressibility  equations  are 
£,  Eq.  (3.17),  and  (pT'),  Eq.  (3.13).  Since  the  equation  of 
continuity  is  in  the  form 

+  =  o  , 

9x  dz  ’ 

Crowley's  second  order  scheme  for  advection  was  modified  to 
use  (pu)  and  (pw)  as  pseudo-velocities. 

This  code  is  now  complete  and  problems  will  be  run  in 
the  next  six  months  in  order  to  compare  with  HAIFA  results. 

In  addition,  a  problem  will  be  calculated  using  a  completely 
compressible  code . 
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3.2  MOISTURE 

3.2.1  Derivation  of  Equations 

In  this  section  the  effects  of  moisture  on  the  equa¬ 
tions  for  Boussinesq  fluid  flow  are  discussed.  Frequently, 
atmospheric  water  in  the  form  of  water  vapor,  cloud  water, 
and  precipitation  will  have  important  effects  on  the  charac¬ 
teristics  of  gravity  waves  caused  by  mountains.  Lee  waves 
are  frequently  accompanied  by  clouds  which  can  be  expected  to 
modify  the  stability  of  the  air  through  the  presence  of  the 
latent  heat  of  condensation  which  the  cloud  water  adds  to 
the  air.  In  addition,  there  are  effects  discussed  by 
Orville^*^  of  up-slope  winds  due  to  high-level  heating 
and  evaporation,  but  these  are  not  of  primary  interest  in 
our  investigation.  Consequently,  the  terms  resulting  in 
changes  of  stability  of  the  air  in  which  clouds  are  forming 
are  of  primary  interest. 

Radiative  heating  and  cooling  of  the  air  has  not  been 
taken  into  account  in  this  discussion,  even  though  the  bound¬ 
ary  condition  on  moisture  is  affected  by  it.  Boundary  layer 
effects  at  present  are  largely  omitted  from  HAIFA  but  are 
considered  in  Appendix  E,  It  will  be  beneficial  to  incorpo¬ 
rate  them,  together  with  radiative  terms,  at  a  later  stage  of 
the  code  development  work. 

The  HAIFA  equations  are  to  be  modified  to  include  the 
effects  of  moisture  by  incorporating  the  following  major 
changes : 

(1)  the  momentum  equation  is  modified  by 
including  the  effects  of  moisture  in 
the  buoyancy  term, 

(2)  the  equation  of  state  for  air  is 
changed  to  include  moisture, 
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(3)  the  energy  equation  is  modified  to  in¬ 
clude  energy  changes  equivalent  to  the 
latent  heat  of  water  being  given  to  or 
taken  from  the  air, 

(4)  a  new  equation  is  added  to  account  for 
the  conservation  of  moisture  in  the  air 
excluding  rain  water,  and 

(5)  a  conservation  equation  is  included 
which  expresses  the  rainwater  content 
in  the  atmosphere  including  sources  and 
sinks  at  the  boundaries. 

An  outline  of  the  derivation  of  the  equations  is  given  below. 

The  momentum  equations  and  the  equation  of  state  for 
a  system  with  moisture  are 


du  1  3p 

Ht  '  p 0  ’ 

(3.18) 

II 

•SlS 

-i-2£  +  £-g(i  +  i  +  i  )  , 

p  3z  p  c  XJ 

0  0 

(3.19) 

p  -  pRT(l  +  Er)  . 

(3.20) 

These  equations  have  been  derived  by  Orville,  and 
Ogura  and  Phillips  among  others.  The  energy  equation  is 
now  derived  from  the  first  law  of  thermodynamics  in  order 
to  redefine  some  terms  used  previously,  i.e., 

Tds  -  dh  -  -2-  .  (3.21) 

P 
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For  an  adiabatic  system,  this  can  be  written  as 


dh  =  1  dp  =  1/ap  + 
eft  “  p  3t  p \9t 


(3.22) 


In  the  Boussinesq  approximation,  the  first  two  pressure  terms 
in  brackets  are  zero  and  the  third  term  is  equivalent  to 
-wgp  (assuming  the  hydrostatic  approximation), 


cTt 


-wg  . 


(3.23) 


The  enthalpy  changes  will  include  energy  changes  due 
to  both  advection  of  the  temperature  and  latent  heat  being 
released  or  absorbed  as  the  moisture  in  the  air  changes  phase. 
The  energy  equation  can  thus  be  expressed  as 


3T  .  ,,3T  .  3T 
5t  *  “57  *  w57 


-wT 


where 


T  =  T  + 


Lr 


7=7  +  T' 

o 

T  * 


where  T  is 


Expressing  the  temperature  as 
a  mean  value  which  is  constant  and  T*  represents  all  vari¬ 
ations  of  the  temperature  from  this  space  averaged  quantity, 
i .  e . , 


**  It 

T  =  T  +  x1  ♦  ~ 

°  p 


7  +  7" 

o 


(3.24) 
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£q.  (3.23)  becomes 


3T" 

at 


+ 


-IF -If 


-wr 


(3.25) 


A  diffusion  term  may  be  included  in  the  above  equation  as 
+ktV2TM  where  kt  is  a  constant  coefficient. 

The  vorticity  equation,  derived  from  the  momentum 
and  continuity  equations,  is 

It  +  ul£  +  wlz  *  g(1  +  lc  *  V  Ix^o5 


+ 


(3.26) 


Making  the  Boussinesq  approximation  and  the  further  restric¬ 
tion  that  T  / (Tq  ♦  T')  ■  1  ,  the  equation  is 


in 

at 


+ 


+  *  t_) 

c  r 


3T,. 

W~ 


♦  fr  r  t1  ♦  *c  ♦  *r>  H  ♦ 

o  p 


at. 


3A. 


&  TT  +  Jx~ 


(3.27) 


Eqs.  (3.25)  and  (3.27)  and  the  Poisson  equation  relating  the 
stream  function  and  the  vorticity  replace  Eqs.  (2.7)  and 
(2.10)  in  the  basic  HAIFA  scheme. 

The  equation  of  water  vapor  conservation  is  obtained 
by  deriving  equations  of  total  water  conservation  and  rain 
water  conservation  and  taking  the  difference  between  them. 
The  total  water  conservation  equation  may  be  obtained  by 
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equating  the  sum  of  the  time  derivative  of  the  total  water 
plus  the  diffusion  of  the  water  carried  as  cloud  water  and 
moisture  to  zero.  Mathematically,  this  is  expressed  as 

fjlPQ)  =  -V-(rrtf)  -  V-(p*ctf) 

-  V-  pf  ( tf-VD)A°  dD 


♦  kV2  p  (r+i-c) 


(3.28) 


The  integral  term  on  the  right  side  of  the  equation  repre¬ 
sents  rain  water  advection  and  fallout  as  a  function  of 
droplet  diameter.  Several  authors  including  Orville,  Kess¬ 
ler,  ^  Srivastava, ^  and  Armason,  et,al.v1^  have  derived 
expressions  for  water  droplet  formation  and  precipitation  in 
the  atmosphere.  At  this  stage  in  the  development  of  HAIFA 
with  moisture,  we  have  elected  to  program  the  conservation 
equations  as  derived  by  Orville  with  one  or  two  exceptions 
noted  below  and  will  modify  these  equations  as  we  derive  or 
discover  better  prescriptions  for  each  of  the  terms. 

The  final  equation  for  total  water  conservation  may 
be  expressed  as 


-VV*Q  +  l 


3V, 
r  9  z 


9 1_  . 

+  v  +  -V2 

t  3z  p 


p(*c+r) 


+  9,  V  +  —  — 
r  t  p  9z 


1  9p 


9 1 


+  b  v  ±  +  V  _ _ 

rVt  p  lx  Vt  9x 


(3.29) 
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The  last  two  terms  on  the  right  side  of  Eq.  (3.29)  are  ignored 
by  Orville.  At  the  present  stage  of  our  analysis,  the  magni¬ 
tude  of  these  terms  with  respect  to  others  in  the  equation  are 
is  unknown;  further  study  will  be  made  to  justify  retaining 
them. 

The  final  equation  required  to  complete  our  analysis 
expresses  conservation  of  rain  water  in  the  atmosphere. 

The  change  of  the  rain  water  with  respect  to  time  is  equal  to 
the  advrection  and  fallout  of  the  droplets  plus  a  source  term 
which  expresses  the  conversion  of  cloud  water  into  rainwater, 
the  growth  of  rain  drops  through  coalescence,  and  the  evapo¬ 
ration  of  rain  falling  through  unsatuiated  air.  The  production 
terms  also  have  been  derived  by  the  authors  already  noted.  The 
most  satisfying  expression  seems  to  be  that  derived  by  Orville 
or  Arnason.  For  consistency,  our  original  equations  for  the 
production  term  will  be  equivalent  to  those  arrived  at  by  Or¬ 
ville.  Modifications  will  be  made  where  we  obtain  an  im¬ 
proved  description  of  the  processes  being  undergone  by  the 
water.  The  equation  may  be  expressed  as 


dl-r  dVt  3V,  £  V, 

_ r  _  9  t  _  r  t  3p  _  .  t  _  r  t  op 

St  ’  lr  3x  "  p'~  3x  '  4r  3z  "  p  3z 


3£. 


-  u 


sir  + 


3£. 


oZ 


3£_ 


3£_ 


t  Jz  vt  S3T 


(3.30) 


This  equation,  Eq.  (3.30),  includes  variations  of  p,  £r,  and 
Vt  with  respect  to  tne  horizontal  direction.  Orville  has  ig¬ 
nored  these,  but  for  completeness  and  until  we  can  substantiate 
that  they  are  negligible,  they  will  be  carried  in  our  studies. 
The  production  term,  Pr  ,  is  described  in  detail  in  Appendix  B. 
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The  equation  for  conservation  of  water  vapor  is  found 
by  subtracting  Eq.  (3.30)  from  Eq.  (3.29).  The  result  is 

3V 

|f  ♦  ?-Vq  -  Icq  V’q  -  »r  -  Pr  .  (3.31) 

Eqs.  (2.8),  (3.25),  (3.27),  (3.29),  and  (3.31)  constitute 
the  complete  set  of  equations  to  be  solved  in  HAIFA  with 
moisture . 

3.2.2  Difference  Equations 

The  difference  equations  used  in  HAIFA  with  moisture 
are  formed  in  an  identical  manner  as  those  in  the  basic  HAIFA. 

All  moisture  terms  are  cell  centered  quantities.  The  time 
differences  are  taken  in  the  forward  direction,  the  advection 
terms  are  treated  by  Crowley’s  schemes  and  all  other  terms  are 
centered  in  space  through  appropriate  averaging.  Since  this 
version  of  HAIFA  is  not  thoroughly  checked  out  at  the  time  of 
this  report,  the  finite  difference  equations  will  not  be  pre¬ 
sented  here.  A  complete  listing  of  these  equations  will  be  a 
part  of  the  annual  report  under  this  contract. 

3.3  VARIABLE  ZONING  IN  VERTICAL  DIRECTION 

The  modifications  to  the  basic  HAIFA  code  that  will 
enable  it  to  operate  with  a  mesh  of  variable  spacing  in  the 
vertical  direction  are  examined  in  this  section.  This  modi¬ 
fication  affords  the  ability  to  resolve  more  finely  certain 
areas  without  excessively  slowing  the  computation  by  re¬ 
quiring  fine  zoning  throughout  the  giid.  Modifications  to 
two  routines  of  the  code  are  necessary.  They  are  the  Poisson 
equation  solver,  and  the  vertical  advection  subroutine.  Each 
modification  is  discussed  below. 
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3.3.1  The  Poisson  Solver 

The  use  of  the  Past  Fourier  Transform  in  the  horizontal 
x-direction  imposes  the  limitation  that  the  spatial  interval, 

Ax  ,  be  constant.  In  the  vertical  direction,  however,  the 
solution  of  the  Poisson  equation  is  obtained  by  Gaussian  elim¬ 
ination  and  is  not  limited  to  a  constant  spatial  interval. 

The  Gaussian  elimination  subroutine  of  POISPK  solves 
a  system  of  difference  equations  approximating 


a  $ 


Q 


(3.32) 


The  solution  of  these  equations  is  briefly  outlined  below: 


The  finite  difference  form  of  Eq.  (3.32) 
may  be  written  as  a  tridiagonal  system 

Ai  ^i+1  +  Bi  *  ci  V'i-i  “  Di  • 


(3.33) 


Lotting 


*i 


(3.34) 


which  implies 

=  Ei_i  i(/.  +  ,  (3.35) 

and  substituting  into  the  tridiagonal  sys¬ 
tem,  the  coefficients  E^  and  G^  may  be 
expressed  as 
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IT  ""+  C.  E 


i  i- 1 


(3.36) 


„  _  Di  '  Ci  Gi-1 
1  Bi  ♦  ci  Ei-i 


(3.37) 


The  finite  difference  form  of  Eq.  (3.32) 
for  constant  vertical  zoning  is 


^i+l  ‘  t2  +  «CAz)2]t|/i  ♦  ^»i_1 


(Az)2 


Qi  • 


(3.38) 


and  the  coefficients  ,  Bi  »  ,  an<* 

are  thus  equivalent  to 

At  -  l/(Az)*  , 


B. 


i  -  -o  -  2/(Az)2  , 


C4  •  1/ (Az) 2  , 

Di  •  q.  . 

Using  these  coefficients,  and  can 

be  calculated  and  thus  the  ^  may  be 
solved  for  recursively. 


(3.39) 


With  variable  zoning  the  finite  difference  form  of 
Eq.  (3.32)  becomes 
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♦i*l  '  *i  '  *i  ' 


TzT 


'i  -1 


Azi  +  Azi-1 


-  a'I'i  ■  Qa  , 


(3.40) 


where  the  location  of  4*  and  Az  are  shown  below, 


The  coefficients 
to 


,  C.  ,  and  D. 

*  1/ (  zi(dzi  +  Azi_1)/2) 
--2/ (A*.  Azi_1)  -  a 


are  now  equivalent 


(3.41) 

Cj.  -  l/(Azi_1  (Aza  +  Azi_1)/2) 

Di  "  Qi  ‘ 


The  values  of  and  are  computed  using  the  above 

coefficients  and  is  computed  in  the  same  manner  as  indi¬ 

cated  above  by  Eq.  (3.35), 
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3.3.2  Vertical  Advection 

The  advection  schemes  discussed  previously  are  valid 
for  uniform  zones  only.  The  equivalent  scheme  for  variable 
size  zones  is  derived  below  for  the  Crowley  second  order 
scheme.  It  has  been  incorporated  into  a  version  of  HAIFA 
which  is  currently  being  tested.  The  fourth  order  scheme  will 
be  considered  at  a  later  date. 

The  one-dimensional  advection  equation  in  conservation 
form  may  be  written  for  flow  in  the  z-direction  as  follows: 

It  +  "  0  *  (3.42) 


where  ^  is  a  variable  representing  the  quantity  to  be  ad- 
vected.  Only  the  one-dimensional  equation  need  be  considered 
due  to  the  splitting  technique  used  in  HAIFA. 

In  finite  difference  form,  Eq.  (3.42)  is 


♦ 


n+1 


♦j  -  SIjMj*!  '  CV+)j] 


(3.43) 


The  term  A(v<t>)  .  requires  the  flux  across  the  boundary  of 
the  j  cell  (see  figure  below) . 
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The  flux  at  the  j+1  boundary  may  be  expressed  as  v^+1 
where  ^  represents  the  value  of  the  variable  $  at  that 
boundary.  Assuming  $  to  vary  linearly  between  zone 
centers,  this  flux  may  be!  expressed  as;  -Tj  •'  ‘  ' 

,  y  w  VM**  •TV  V  •'V*  f  >*  *•? ' 

•  .  *n+l  '  V-%  .  '  •  ;• 


boundary  "It  J  n  v*  « "  2 
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Assuming  $  ■  a  ♦  bz  ,  and  integrating, 


(3.44) 
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(3.45) 


The  new  value  of 


can  then  be  expressed  using  Eq.  (3.45) 
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4.  TEST  PROBLEMS 

Several  problems  have  been  calculated  using  the  basic 
HAIFA  code.  The  results  of  each  are  presented  in  this  sec¬ 
tion  and  comparisons  with  other  results  are  made  where  pos¬ 
sible.  An  edit  routine  to  determine  the  momentum  flux 
(wave  drag  associated  with  gravity  waves)  was  written  and 
is  described  in  detail  in  Appendix  C. 

Table  I  summarizes  the  initial  conditions  used  for 
each  problem.  The  boundary  conditions  in  each  case  were 
those  described  in  Section  2.4  of  this  report.  The  grid 
size  consisted  of  35  vertical  cells  by  64  horizontal  cells. 

•  .  *  / 

4.1  SINGLE  WAVE 

The  atmospheric  and  horizontal  velocity  conditions 
to  produce  a  single  gravity  wave  were  arrived  at  using  the 
results  presented  on  two-dimensional  mountain  lee  waves  by 
Palm  and  Foldvik.^11^  They  had  established  that  if  the 
quantity 

S  1  3*u  - 

u*’  ~  u  3z2  ' 

where  S  is  the  stability  of  the  atmosphere,  has  a  value 
at  the  ground  level  which  is  at  least  2.5  times  as  large  as 
the  minimum  value  (usually  located  7-10  km  above  the  ground), 
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Figure  4.2  -  Two  Wave  Velocity  Profile. 
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the  wave  motion  in  the  lower  troposphere  depends  only  on  the 
wind  profile  and  the  stability.  This  condition  is  almost 
always  satisfied  when  mountain  waves  occur.  A  diagram  giv¬ 
ing  the  expected  wave  lengths  of  lee  waves  under  various 
stability  and  wind  profiles  was  presented.  In  particular, 
regions  of  one  and  two  waves  were  indicated.  Using  this 
diagram,  a  single  wave  of  approximately  16  km  in  length  was 
predicted  for  a  lapse  rate  equal  to  one-hrlf  the  dry  adiabatic 
value  (see  Figure  4.3),  and  the  exponential  velocity  profile 
shown  in  Figure  4.1. 

The  numerical  results  calculated  using  HAIFA  are 
shown  in  Figures  4.4  through  4.7  as  streamlines  and  vertical 
velocity  contours  at  several  times  up  to  1-1/4  hours.  The 
measured  wave  length  from  Figure  4.5  or  4.7  is  approximately 
15  km.  As  can  be  observed  from  the  results,  only  one  wave 
did  form  during  the  time  the  problem  was  run.  The  cyclic 
boundary  condition  prevented  any  further  computation  due  to 
disturbances  created  by  the  obstacle  in  the  flow  stream 
being  introduced  into  the  main  flow  upstream  of  the  mountain. 
Some  interference  with  the  upper  boundary  positioned  at 
10.9  km  may  also  be  seen  at  the  latest  times. 

The  momentum  edits  u’v*  (see  Appendix  C)  located 
one  cell  or  312.5  meters  above  the  mountain  top  are  shown 
in  Figures  4.8  and  4.9  for  various  lengths  used  in  obtaining 
the  horizontal  averages.  The  qualitative  result  obtained 
from  these  figures  indicates  a  decrease  in  the  edited  quan¬ 
tity  as  the  length  used  in  the  averaging  length  is  increased, 
i.e.,  a  lower  amount  of  drag  is  created  by  the  mountain.  One 
exception  appears,  however;  this  can  be  noted  as  a  cross  over 
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of  two  of  the  curves  occurring  at  approximately  3400  to  4000 
seconds  on  either  of  the  figures.  The  same  phenomenon  occurs 
when  the  averaging  length  is  reduced  by  discounting  zones 
from  in  front  of  the  obstacle  as  well  as  the  rear.  While 
it  is  not  clear  what  the  averaging  length  should  be  in  these 
cases  or  the  intrepretation  of  these  results ,  it  is  clear 
that  the  magnitude  of  the  edited  quantity  is  only  equal  to 
the  drag  on  the  mountain  if  the  inlet  and  outlet  values  of 
p  and  pu2  are  identical.  Since  this  is  the  case  only 
when  the  total  numerical  grid  length  is  used  as  the  averaging 
length,  due  to  the  cyclic  boundary  conditions,  a  value  for 
the  drag  on  the  mountain  can  only  be  estimated  from  the 
uppermost  curve  of  Figure  4.8.  The  value  of  the  drag 
reached  at  4,44S  seconds  was  approximately  equal  to  10  dynes/ 
cm2.  This  value  agrees  qualitatively  with  measured  values 
of  the  momentum  flux  reported  by  D.  K.  Lilly^12^  for  mea¬ 
surements  at  Boulder,  Colorado. 

One  other  important  feature  of  the  momentum  flux 
edit  is  the  oscillatory  character  of  the  values  with  time. 
This  is  thought  to  be  related  to  the  formation  of  the  in¬ 
dividual  vertical  velocity  cells,  i.e.,  as  a  new  positive 
or  negative  cell  is  formed,  the  effect  seems  to  be  to  in¬ 
crease  or  decrease  the  horizontal  average  of  the  vertical 
flux  of  horizontal  momentum.  This  cyclic  character  is 
perhaps  more  clearly  seen  in  the  edits  of  the  two  wave 
problem  discussed  later. 

Figure  4.10  shows  the  momentum  edit  as  a  function 
of  height  at  a  time  of  4,445  seconds.  The  value  goes  to 
zero  very  quickly  above  the  mountain.  This  indicates  the 
solution  is  not  yet  approaching  a  steady  state  value  since 
the  drag  for  a  steady  problem  would  be  constant  with  height. 
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4.2  TWO  WAVE  PROBLEM 

The  conditions  for  the  two  wave  problem  closely  match 
those  of  a  problem  described  by  Wurtele  and  Foldvik. 

These  authors  computed  numerically  the  transient  formation  of 
a  mountain  lee  wave.  One  wave  of  length  10  to  15  km  was 
formed  using  conditions  similar  to  those  described  for  the 
two  wave  problem  in  Table  I.  Using  the  same  type  of  analysis 
as  described  previously  for  the  single  wave  problem,  the  Palm 
and  Foldvik  results  indicated  a  wave  of  approximately  9.2  km 
wavelength  as  well  as  a  second  wave  of  approximately  25  km 
wavelength  should  be  present.  Private  communication  with 
Wurtele  indicated  that  the  longer  wave  was  not  noticed  in 
their  calculation.  These  data  were  incorporated  into  HAIFA 
and  run  to  a  time  of  5,474  seconds.  The  results  of  this 
numerical  calculation  are  presented  in  Figure  4.11  through 
4.14.  The  streamlines  shown  in  Figure  4.12  at  a  time  of 
5,474  seconds  display  the  presence  of  two  waves.  The  shorter 
wave  appears  just  above  and  behind  the  obstacle  and  agrees 
with  theory  in  that  the  wave  length  is  approximately  12  km  in 
length.  A  second  wave  appears  behind  the  obstacle  at  a 
height  of  7  to  8  km  with  a  wavelength  of  approximately  22  km. 
The  sma’.l  discrepancies  between  predicted  and  calculated 
wavele.  gths  are  probably  the  result  of  the  linear  theory 
used  in  producing  the  diagram  of  Palm  and  Foldvik.  The 
presence  of  the  second  wave  is  also  strongly  evident  in  the 
bending  over  of  the  vertical  velocity  fields.  The  interaction 
of  the  two  waves  is  seen  by  the  presence  of  small  vertical 
velocity  regions  whose  direction  is  opposite  of  the  velocity 
cell  completely  surrounding  it. 

The  figures  showing  the  growth  of  the  vertical  veloc¬ 
ity  cells  may  be  compared  to  the  results  of  Wurtele  and  Fold¬ 
vik.  The  Figure  4.15  was  taken  from  their  paper  and  indicates 
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Figure  4. IS  -  The  field  of  vertical  velocity  as  cal¬ 
culated  by  Wurtele  and  Foldvik.  The 
lowest  panel  represents  the  total 
streamline  field. 
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a  large  negative  vertical  velocity  cell  over  the  obstacle  for 
all  time.  The  positive  cell  just  downstream  of  the  obstacle 
is  moving  toward  the  upstream  side  of  the  obstacle  with  time 
but  does  not  ever  dominate  the  flow  over  the  obstacle.  Our 
calculations  indicate  a  negative  cell  over  the  obstacle  for 
times  to  S,300  seconds.  Shortly  after  that  time,  the  large 
positive  cell  at  the  rear  of  the  obstacle  combines  with  the 
small  positive  cell  forward  of  the  obstacle.  A  small  nega- 
tive  cell  still  remains  over  the  obstacle  itself,  but  there 
is  no  longer  a  dominate  negative  or  downwind  flow  over  the 
obstacle.  These  differences  are  not  well  understood.  The 
phenomena  may  be  due  to  differences  in  initial  or  boundary 
conditions  as  none  of  these  were  completely  indicated  in  the 
Wurtele  and  Foldvik  paper.  The  additional  time  of  the  com- 
putation  or  the  presence  of  the  second  wave  may  also  be 
responsible  for  this  phenomena.  The  length  of  the  obstacle, 
which  was  not  specified  in  the  previous  study,  was  found  in 
HAIFA  calculations  to  have  a  definite  influence  on  the  pat¬ 
tern  of  vertical  velocities  over  the  obstacle.  The  tempera¬ 
ture  distribution,  which  was  specified  as  equal  to  one-half 
the  dry  adiabatic,  may  also  vary  a  small  amount.  Any  of 
these  parameters  may  have  caused  the  small  differences  seen 
in  our  calculation  and  that  reported  by  Wurtele  and  Foldvik. 
The  similarities,  particularly  of  the  vertical  velocity  cells 
at  times  less  than  5,000  seconds,  certainly  indicate  good 
qualitative  agreement. 

The  momentum  edit  of  the  two  wave  problem,  shown  in 
Figure  4.16,  indicates  a  much  more  cyclic  nature  than  of  the 
single  wave  problem.  Each  half  cycle  appears  to  have  some 
correspondence  to  the  complete  formation  of  a  vertical  cell 
although  the  interaction  of  the  two  waves  present  makes  the 
intrepretation  of  the  data  difficult. 


73 


3SR-79S 


4.3  UNIFORM  VELOCITY 

A  problem  using  a  velocity  distribution  uniform  with 
height  and  equal  to  10  m/sec  perturbed  by  a  one  kilometer  high 
mountain  was  completed.  The  lapse  rate  was  set  equal  to  one* 
half  the  dry  adiabatic.  Figures  4.17  through  4.20  show  the 
resulting  streamlines  and  vertical  velocity  cells  formed  under 
these  conditions  Wurtele  and  Foldvik  have  also  invest!* 
gated  this  problem,  the  results  of  which  are  shown  in 
Figures  4.21  through  4.23.  A  comparison  of  their  streamlines 
with  our  results  show  a  continuous  spectrum  of  waves  is  ex¬ 
cited  in  both  calculations,  which  when  added  together  produce 
growing  numbers  of  upwind-tilting  troughs  and  crests  ex¬ 
tending  to  great  heights.  The  figures  showing  the  verti¬ 
cal  velocity  cells  at  the  forward  and  rear  of  the  obstacle 
show  these  upwind-tilting  troughs  and  crests  even  more  dis¬ 
tinctly.  Lyra^1*^  theoretically  showed  these  same  results 
using  a  linear  analysis.  His  steady  state  analytical  result 
for  the  streamlines  and  the  vertical  velocity  field  are 
shown  in  Figures  4.24  and  4.2E.  While  there  are  certainly 
similarities  in  the  results  of  Lyra,  Wurtele  and  Foldvik,  and 
the  S1  calculation,  there  are  also  some  significant  differ¬ 
ences.  The  four  total  streamline  fields  computed  by  Wurtele 
and  shown  in  Figures  4.21  and  4.22  show  a  large  amplitude 
wave  just  above  the  lee  slope.  The  vertical  velocity  in  this 
region  is  more  than  five  tiroes  the  upstream  wind  and  the 
total  horizontal  velocity  is  negative  at  some  grid  points. 

This  feature  is  not  present  in  the  linear  theory  and  did  not 
appear  in  the  S1  computations.  We  are  continuing  to  investi¬ 
gate  these  differences. 

One  of  the  most  significant  items  found  in  calculating 
this  problem  was  a  numerical  instability  associated  with  the 
flow  when  the  normal  stability  criteria  for  the  advective 
terms  of  the  equations  was  used.  An  initial  computation 
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Figure  4,19  -  Vertical  Velocity  Field  from 
Uniform  Velocity  Problem. 
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Figure  4.23  -  Field  of  vertical  motion  computed  under 
upwind  conditions  for  successive  times 
from  Uniform  velocity  problem. 


Figure  4.24  -  Streanlines  fro*  the  linear  theory 
(after  Lyra). 


Figure4.25  -  Field  of  vertical  velocity  when 
U  -  constant  with  height  (after 
Lyra).  Isopleths  for  w  >  0  only. 
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using  this  tine  step  control  produced  a  series  of  large  wave 
length  high  amplitude  waves  which  propagated  throughout  the 
flow  very  quickly.  The  problem  was  recalculated  by  putting 
an  upper  limit  on  the  time  step  which  was  based  on  the  phase 
speed  of  the  largest  of  these  waves,  i.e.,  a  wave  with  a 
50  km  wavelength.  This  limited  the  time  step  to  less  than 
14  seconds  per  cycle  assuming  the  overall  criteria  to  be  that 
the  SO  km  wave  would  not  completely  traverse  a  grid  cell  in 
one  cycle.  The  actual  limiting  time  step  used  in  the  recal- 
culation  was  12.0  seconds.  The  resulting  wave  pattern  is 
the  one  shown  and  previously  discussed  in  this  section. 

This  new  stability  criteria,  which  had  not  been  previously 
used,  was  not  required  in  earlier  problems  due  to  either 
(1)  the  damping  of  the  disturbances  caused  by  the  wind  shear 
or  (2)  the  high  velocities  in  the  single  wave  problem  con¬ 
trolling  the  time  step  to  an  acceptable  value.  Later  prob¬ 
lems,  the  tropopause  and  inversion  layers,  exhibited  this 
same  instability. 

4.4  INVERSION  LAYER  I 

The  determination  of  the  effect  of  an  inversion  layer 
in  the  atmosphere  was  calculated  using  the  basic  HAIFA  -.ode. 
The  inversion  layer  was  described  as  a  positive  4#C  tempera¬ 
ture  change  over  a  1.5  km  height  as  shown  in  Figure  4.26. 

The  other  initial  conditions  are  described  in  Table  I.  The 
results,  shown  in  Figures  4.27  through  4.30,  indicate  a 
small  effect  in  the  vertical  velocity  cells  at  heights  cor¬ 
responding  to  the  inversion  heights.  The  cells  appear  to  be 
broader  at  a  5  km  height  than  those  seen  in  the  two  wave 
case  for  example.  There  also  appear  to  be  displacements  in 
the  vertical  cells  at  this  position.  However,  these  may  be 
due  more  to  the  change  in  the  lapse  rates  at  this  position 
than  the  presence  of  the  4°C  temperature  increase. 
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Figure  4.26  -  Temperature  Distribution  used  in  Inver> 
sion  Layer  Problem. 
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Figure  4.29  -  Vertical 
version 
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Because  of  the  coarse  zoning  at  the  inversion  layer, 
the  definition  of  the  flow  is  poor.  This  problem  will  be 
recalculated  using  the  variable  zoning  code  and  the  results 
will  be  reported  in  the  final  report  of  the  contract. 

4 . S  TROPOPAUSE  PROBLEM 

The  test  calculation  representing  a  tropopause  prob¬ 
lem  consisted  of  initial  conditions  as  described  in  Table  I. 
The  calculated  streamlines  and  vertical  velocity  contours 
are  shown  in  Figures  4.31  through  4.34,  The  most  noticeable 
characteristic  of  the  resulting  solution  is  the  tilting  of 
the  vertical  velocity  cells  toward  the  upwind  direction. 

The  streamline  pattern  for  this  problem  did  indicate,  but 
not  clearly,  this  same  phenomena  of  the  upwind  tilting  of 
the  gravity  wave  peaks. 
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5.  RADIATION  IN  THE  EARTH'S  ATMOSPHERE 

The  radiative  transfer  problem  in  the  Earth's  atmos¬ 
phere  reduces  to  the  solution  of  the  seemingly  simple  equa¬ 
tion 


dIv.  j 
3s”  Jv 


(5.1) 


which  states  that  radiant  intensity,  in  traversing  the  element 
of  length  ds  ,  will  be  augmented  by  sources  in  the  amount 
Jv  ds  and  diminished  by  extinction  in  the  amount  icvIvds  • 

In  general,  Iv  ,  the  radiant  intensity,  and  Jy  ,  the  source 
function,  depend  on  both  a  spatial  coordinate  r  and  an  angular 
coordinate  (direction)  ft  at  the  point  r  ,  as  well  as  upon 
the  frequency  v  .  The  time  dependence  of  these  quantities  is 
ignored  because  the  radiative  state  of  the  atmosphere  is 
established,  for  all  practical  purposes,  instantaneously.  <v 
is  the  extinction  coefficient,  which  describes  the  relative 
depletion  in  the  intensity  of  the  beam,  d I v/ I v  ,  upon  tra¬ 
versing  the  element  of  distance  ds .  <v  is  in  general  the 
sum  of  an  absorption  part  and  a  scattering  part.  Jv  describes 
the  additions  made  to  the  beam  intensity  along  ds  by  thermal 
or  non-thermal  emission  and  by  scattering.  In  the  case  where 
the  source  consists  only  of  thermal  emission,  Jy  does  not 
depend  on  Iv  and  Eq.  (S.l)  can  be  solved  explicitly  for  Iv: 
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Iv(s,8) 


exp 


%(*’)  ds' 


(5.2) 


♦ 


fs 

r  rs 

1  J  (s",8)  exp 

-  j  *v(s')  ds' 

\ 

;S" 

k 

ds" 


where  sq  corresponds  to  some  boundary  at  which  1^  is  pre¬ 
sumed  known.  However,  even  in  the  case  where  Jv  depends  on 
Iv  ,  Eq.  (S.2)  is  a  perfectly  valid  alternate  formulation  of 
the  radiative  transfer  equation.  Eq.  (S.l)  will  be  called  the 
differential  form,  and  Eq.  (S.2),  the  integral  form,  of  the 
transfer  equation. 

The  consideration  will  be  limited  to  plane  geometry, 
which  is  justified  in  view  of  the  small  depth  of  the  Earth's 
atmosphere  in  comparison  to  its  radius  (when  the  sun  is  near 
the  horizon  the  plane-parallel  atmosphere  approximation  fails). 
The  geometry  is  illustrated  in  Figure  5.1(a).  The  vertical 
coordinate  is  denoted  by  z  ,  and  the  angular  coordinates  de¬ 
fining  the  direction  ft  at  z  are  denoted  by  6  and  $  . 

As  is  customary,  the  variable  m  =  cos8  is  employed  rather 
than  8  itself.  From  Figure  S.l(b),  it  is  clear  that  in  these 
coordinates  the  element  of  distance  ds  is  related  to  dz  by 


so  that  Eq.  (S.l)  becomes 

3Iv 

»TT  "  Jv  "  S  *v 


(5.3) 
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*  (a)  (b) 

Figure  5.1  -  Coordinate  system  for  radiation  problem. 

where  Iv  ■  Iv(z,u,4)  and  •  kv(z)  .  If  we  sake  the 
assumption  of  local  thermodynamic  equilibrium  (LTE) ,  which 
is  valid  below  about  70  km  in  the  atmosphere,  then  the  source 
term  Jy  may  be  replaced  by  a  more  explicit  expression, 
leading  to  the  transfer  equation 

4r  ‘  «J»v  •  V  *  6v  77/ V*-8-8')  V*-8'1  da’  -  Jv  (5-4> 

►  * 

where  o^  is  the  volume  absorption  coefficient  ov  corrected 
for  stimulated  emission 

*  «v(l  -  e’hv/kT)  ,  (S.5) 

Pv  is  the  scattering  coefficient,  Bv  is  the  Planck  function 

VT>  •  ^7/^77  •  <5-6> 

and  Pv  is  the  phase  function,  defined  so  that 

t>vU,8,8')  $ 
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is  the  probability  that  a  photon  entering  a  volume  element 
around  z  from  direction  ft'  will  be  scattered  into  the 
cone  dfl  of  directions  around  ft.  S  nee  absorption  is  ex* 
plicitly  represented  in  Eq.  (5.4),  the  above  probabilities 


must  sum  to  one 


/p  C*.a.S)  - 

Jkm 


(5.7) 


rather  than  to  some  number  less  than  one  as  would  be  the  case 
if  absorption  were  tacitly  included  in  the  scattering  terms. 

Because  a  volume  element  in  the  atmosphere  has  no  pre¬ 
ferred  direction  with  respect  to  scattering  (due,  say,  to  a 
permanent  dipole  moment),  the  scattering  probability  function 
Pv  depends  only  on  the  angle  between  ft  and  ft'  .  If 


"s  •  eoses 


this  means  that  P  ■  PU(*,WC)  •  is  expressed  in  terms  of 

known  angles  6  ,  8*  ,  4  ,  V  as  follows: 

ws  ■  ft*ft'  ■  (sin8  cost,  sine  sint,  cos8) 


•  (sine*  cost',  sine*  sint',  cos8') 


•  sine  sine' (cost  cost'  ♦  sint  sint') 


♦  cose  cose* 


yy'  ♦  /l  -  y 2/l  -  y*2  cos  (t"t') 


Let  us  integrate  Pv  over  all  azimuthal  directions  t  ; 
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*  /lV1  costjdt 


The  second  and  third  integrals  cancel  one  another,  because  of 

a 

the  periodicity  of  cost  >  leaving 

fU 

pv(z,p,u')  =  J  PvU,u5)dt 

2w 

Pv(z,mm'  ♦  /l-p2  /lV2  cost)dt  .  (5.8) 
o 

The  important  point  to  note  here  is  that  Pv  does  not  depend 
on  4> '  . 

Using  Eq.  (5.8),  one  may  integrate  Eq.  (5.4)  over  azimuth 
t  and  so  deal  only  with  the  azimuthally-averaged  intensity  Iv  : 
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31 


Vt-z-  •  o’  (B  -  I  ) 
h3z  vv  v  v' 


♦  B. 


PvCz»p»p,)iv(2»p,)dy'  -  i, 


(5.9) 


where 


Iv(z.p) 


Iv(z»W»$)d4> 


In  problems  such  as  the  solar  aureole,  the  location  of  the 
neutral  points  in  the  sunlit  sky,  etc.,  it  is  clearly  neces¬ 
sary  to  retain  the  ^-dependence  of  the  intensity;  for  the 
computation  of  some  important  angular  moments  of  the  intensity, 
however,  and  in  particular  for  the  computation  of  vertical 
radiative  fluxes  and  hence  heating  rates,  Iv  contains  all 
the  necessary  information. 

No  more  than  the  first  three  angular  moments  of  the 
intensity  will  be  considered  in  what  follows.  They  are,  in 
the  customary  notation, 


£VU)  *  3  /*Iv(z,fi)dn  m—f  Iv(z,jj)dp 

J  -1 

».,„(*>  ■  A,  Iw(*.8)«  (5.10) 

-kfa a  S  V*-*)" 
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where 

o  .  /l-y2  cos4)  ,  n  »  /l-y2  sinij)  ,  SI  -  y 
a  y  z 

E  ,  F  ,  and  P  may  be  interpreted  physically  as  the  density 
of  radiant  energy,  flux  of  radiant  energy,  and  radiation  pres¬ 
sure  tensor,  respectively.  Two  moments  of  special  interest 
are  the  vertical  flux  Fz  and  the  zz-component  of  pressure 
P 

*zz  ' 


-  fv  lvU,8)dn 


lv(z,y)<*y 


PVU)  ” 


I  y>iv(z,fi)dn 

f  Y  y2^v(z*y)dy 


C5.ll) 


where  for  conveience  the  coordinate  subscripts  have  been 
omitted.  The  other  components  of  the  flux,  F  and  F  , 
might  be  of  interest  for  some  applications,  such  as  the  heat¬ 
ing  of  inclined  slopes,  but  their  calculation  would  require 
retaining  the  ^-dependence  of  Iv  ,  as  would  the  calculation 
of  any  of  the  other  pressure  components. 
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In  general,  radiation  energies  and  pressures  within 
the  atmosphere  are  completely  negligible  compared  tc  material 
energies  and  pressures,  whereas  radiative  fluxes  are  compar¬ 
able  to  other  energy  fluxes  in  the  atmosphere  (latent  heat, 
sensible  heat,  etc.);  the  reason  for  this  i>  that  the  rela¬ 
tively  small  amounts  of  radiant  energy  travel  at  the  speed  of 
light,  while  the  speed  of  material  energy  propagation  is 
essentially  limited  by  the  sound  speed. 

The  definitions  of  radiation  energy  and  flux  given 
above  are  made  plausible  by  looking  at  the  result  of  inte¬ 
grating  Eq.  (5.9)  over  all  y  (remembering  the  normalization 
Eq.  (5.7)  of  Pv): 


dF 

v 

zr 


°i(4*Bv 


cV 


(5.12) 


The  terms  on  the  right-hand  side  are  source  and  sink  terms  to 
the  radiation  energy  field;  if  they  cancel,  then  Eq.  (5.12), 
with  the  above  interpretation  of  Fv,  becomes  the  usual  ex¬ 
pression  for  steady-state  radiative  energy  conservation.  If 
they  do  not  cancel,  then  clearly  more  energy  is  entering  an 
infinitesimal  horizontal  layer  than  is  leaving  it,  or  vice 
versa,  and  the  deposited  or  withdrawn  energy  will  result  in 
a  net  heating  or  cooling  of  that  layer.  The  expression  for 
the  heating  rate  is,  in  fact,  taking  the  origin  of  z-coordinates 
at  the  surface  of  the  Earth 


dT  _  _  dF 
cTt  "  ~  <Tz  * 


(5.13) 


where 
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F(z)  -  /  F  (z)dv  (5.14) 

Jo 

p(z)  ■  density  of  air  at  z 

C  *  specific  heat  of  air  at  constant 
"  pressure  . 

Eq.  (5.13)  is  simply  a  restatement  of  the  first  law  of  thermo¬ 
dynamics,  which  is,  in  the  usual  notation, 

dq  «  de  +  p  dv 

-  Cv  dT  ♦  pRT  d(i) 

-  Cv  dT  -  SI  dp 

for  an  ideal  gas.  Since  any  atmospheric  process  at  a  given 
level  z  will  for  all  practical  purposes  take  place  at  con¬ 
stant  pressure, 


dp  ■  0  ■  R(p  dT  +  T  dp) 
then  the  first  law  can  be  written 

dq  -  Cy  dT  -  S(.p  dT) 

-  (C„  ♦  R)  dT  ■=  C  dT  . 
v  p 

Dividing  by  dt  ,  and  noting  that  the  heating  rate  due  to 
radiation  corresponds  to  ,  leads  to  Eq.  (5.13) 

(the  factor  p  converts  heating  per  unit  mass  to  heating 
per  unit  volume,  in  agreement  with  the  definition  of  F  ). 
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To  obtain  the  radiative  heating  rates  in  Eq .  (5.13) 
is  the  primary  goal  of  our  calculation.  To  arrive  at  these 
numbers,  Eq.  (5.9)  (or  its  corresponding  integral  form)  must 
be  solved  for  the  Iv’s  ,  which  must  be  integrated  ac¬ 
cording  to  Eq.  (5.11)  to  obtain  the  F^'s  ,  and  finally  the 
Fv's  must  be  integrated  in  accordance  with  Eq,  (5.14)  to 
obtain  the  F*s  .  The  mechanics  of  solving  Eq.  (5.9)  will 
be  elaborated  upon  in  the  remainder  of  this  section. 

The  solution  of  Eq.  (5.9)  can  be  separated  into  sev¬ 
eral  sub  ■'tasks,  which  are: 

(1)  specification  of  the  phenomenological 
parameters  entering  the  equation; 

(2)  specification  of  boundary  conditions; 
and 

(3)  discretization  of  the  independent  vari¬ 
ables  y  ,  z  ,  v  for  numerical  solution. 

These  three  sub -tasks  are  discussed  in  turn  in  Sections  5.1 
through  5.3. 

5.1  PARAMETER  SPECIFICATION 

The  parameters  required  in  Eq.  (5.9)  are  the  absorp¬ 
tion  coefficient  ,  the  scattering  coefficient  0V  ,  and 
the  $- averaged  phase  function  Pv  ,  and  6V  are  known 

to  have  both  a  temperature  (T)  and  pressure  (p)  dependence, 
so  that  the  temperature  and  pressure  structure  of  the  atmos¬ 
phere  constitutes  required  input  data  (perhaps  from  a  GCM) . 
Other  pertinent  input  data,  required  for  the  computation  of 
0V  and  Pv  ,  are  the  aerosol  and  cloud  structure  of  the 
atmosphere;  more  specifically,  the  number  density  of  aerosol 
particles  (including  cloud  droplets)  as  a  function  of  both 
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height  and  particle  radius,  and  the  frequency  dependent  index 
of  refraction  of  all  aerosol  constituents.  For  the  computa¬ 
tion  of  optical  path  lengths,  the  mixing  ratios  of  the  non- 
uniformily  mixed  gases  ^0  and  0^  as  a  function  of  height 
are  also  required. 

Needless  to  say,  the  specification  of  atmospheric 
structure  in  such  detail  is  beyond  the  capabilities  of  either 
experiment  (soundings)  or  theory  (GCM's)  at  the  present  time. 
However,  experimentation  with  the  detailed  model  being  devel¬ 
oped  should  point  the  way  toward  simpler  specifications  of 
structure  which  are  nevertheless  sufficient  for  computing 
heating  rates.  At  the  same  time,  one  may  expect  an  improve¬ 
ment  in  vertical  resolution  and  aerosol  prediction  capability 
in  the  GCM's  and  more  accurate  experimental  data,  particularly 
with  regard  to  aerosols,  in  the  near  future.  Hopefully,  this 
will  lead  to  a  felicitous  convergence  of  the  radiation  model's 
need  for  atmospheric  structure  data  and  the  ability  of  theory 
and/or  experimental  to  furnish  it. 

The  absorption  coefficient  a v  will  not  be  discussed 
in  detail  here.  Absorption  in  the  atmosphere  takes  place  in 
large  part  in  vibration-rotation  bands  of  ^0,  CC^,  0^  and 
other  minor  constituents;  each  band  contains  thousands  of 
spectral  lines  ,  resulting  in  an  extremely  rapid  variation  of 
av  with  frequency.  Voluminous  compilations  of  exist, 

but  it  is  impossible,  for  reasons  of  computer  storage  and 
economy,  to  discretize  I  in  frequency  space  finely  enough 
to  follow  the  variations  of  ay  ;  spectral  intervals  must  in¬ 
stead  be  taken  which  contain  many  lines.  Hence,  this  whole 
discussion  falls  more  naturally  into  the  province  of  Section 
5.3,  where  discretization  in  v-space  is  discussed. 

The  Rayleigh  scattering  coefficient  $  ,  D  and  phase 

v  $  **  ns'i 

function  (independent  of  v)  are  discussed  by  Penndorf.  J 
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He  gives  the  index  of  refraction  ng  of  air  at  p  *  760  nun  Hg  , 
T  *  15°C  ,  and  water  vapor  pressure  f  (in  mm  Hg)  as 


(ns  -  1)  x  10*‘  -  64.328  +  29498. l|l46  -  p) 


-  i 


+  255.4 


(41'Xt) 
-  (»••««  - 


where  X  is  the  vacuum  wavelength  in  microns.  The  Rayleigh 
volume  scattering  coefficient  is  then 


8n! 


t"s  - 
»*»; 


IN 


(5.15) 


and  the  Rayleigh  phase  function 


(16) 


is 


PR^ 


W’nJ  /. 
rmr  ^ 


(5.16) 


where 

s 


cos0s  =  cosine  of  scattering  angle 


=  number  density  of  air  molecules  at 
s  760  mm  Hg  and  15°C  -  2.54743  x  1019  cm'9 

Pn  =  depolarization  factor 

N  -  number  density  of  air  molecules  at  p  and 
T  of  interest 


The  factor  involving  pn  expresses  the  effect  of  the  optically 
anisotropic  molecules  upon  the  scattering,  and  its  value  has 
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been  calculated  and  measured  by  a  number  of  investigators 
since  the  effect  was  first  discovered  (see  Penndorf,  Table  II). 
Penndorf  chooses  pn  *  0.035  ,  which  we  shall  use  in  our  work. 
The  number  density  N  may  be  evaluated  from  the  perfect  gar 
law  as 


N  "  fa 

where  k  is  Boltzmann's  constant,  T  is  in  8K,  and  p  is 
in  compatible  units. 

The  scattering  from  aerosols  in  the  atmosphere  (in¬ 
cluding  clouds)  may  be  treated  by  the  Mie  theory.  This  in¬ 
volves  a  certain  degree  of  approximation,  in  that  the  Mie 
theory  assumes  spherical  particles  and  natural  atmospheric 
aerosols  may  not  be  spherical  (although  water  droplets,  the 
most  important  aerosols,  are  indeed  spherical  provided  they 
do  not  have  an  appreciable  fall  velocity).  Also,  the  complex 
indices  of  refraction 


ra  =  n  -  i  n  (5 . 17) 

1  2 

of  the  aerosols,  which  are  required  by  the  Mie  theor/,  are  in 
many  cases  not  well-known,  especially  in  the  infrared 
Nevertheless,  Mie  theory,  despite  its  well-known  mathematical 
and  computational  complexities,  is  the  only  reasonable  approach 
to  aerosol  scattering  currently  available. 

If  a  is  the  radius  of  a  single  spherical  particle, 
having  index  of  refraction  m  relative  to  the  surrounding 
medium,  then  the  scattering  pattern  of  that  particle  for  light 
of  wavelength  X  can  be  described  in  terms  of  the  following 
two  functions : 
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i 

i 


E  HTH?TT(yo'm>Vl,s>  *  bn(“"m>V,s>] 


i 

2 


E  HTKTrr[an(“>m)Tn^s>  *  bn C« ,m3 Tr„ (lJs >] 


(5.18) 


(5.19) 


where 


2ira 

a  "  nr 

(5.20) 

8  ■  ma 

(5.21) 

“s  •  coses 

an  ’  «n{c»)tJC8)  -  a  *n(6)c^ta) 

(5.22) 

♦nc*)*i(«)  -  m  M^n^ 
bn  '  *nlB)cJl«)  -  ■  tn(o)*i(S) 

(5.23) 

"nOO  -  p^Cw) 

(5.24) 

Tn (m)  c  M7Tn(p)  -  (l-y2)n^(y)  . 

(5.25) 

The  PR  are  Legendre  polynomials.  The  ^n»^n  are  called 
Ricatti -Bessel  functions,  and  are  defined  in  terms  of  the  more 
familiar  spherical  Bessel  functions  jn  ,  y  ,  etc.  by 
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*n(2)  "  zjn^z) 

xnfO  ■  -*ynU) 

Cnc^)  ■  *  iXn  ■  zhn2)(z)  • 

For  unpolarized  incident  light,  the  distribution  of 
scattered  intensity  from  the  spherical  particle  is  propor¬ 
tional  to  i  ♦  i  .An  unpolarized  incident  beam  is  also 
12 

tacitly  assumed  in  the  form  (5.16)  of  the  Rayleigh  scatter¬ 
ing  pattern.  The  full  polarization-dependent  treatment  of 
radiative  transfer,  in  which  the  intensity  is  replaced  by  a 
four-component  vector  and  all  the  phase  functions  are  re¬ 
placed  by  4x4  phase  matrices,  involves  a  great  deal  more 
computing  than  the  present  method  for  a  relatively  small  im- 
provement  in  accuracy^  (the  largest  reported  errors  in  Iv 
from  neglecting  polarization  are  about  10  percent,  with  more 
typical  values  being  1-5  percent).  Therefore,  the  present 
model  ’ivill  be  constructed  assuming  every  Mie  scattering  event 
produces  only  unpolarized  light,  of  intensity  proportional  to 
i  +  i 

l  2 

The  extinction,  scattering,  and  absorption  cross- 
sections  for  the  spherical  particle  may  be  computed  to  be 

00 

°ext  ‘  I?  £  (2ntl)  Re(W 
n=l 

00 

°sca  ’jfZ  (2n*l)(|an|2  ♦  |bn|*) 
n=  1 


a  .  =o  ^  -  a 
abs  ext  sea 
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Presuming  that  the  atmospheric  aerosol  at  a  given 
height  contains  a  number  density  N  of  spherical  particles, 

ov  r 

with  a  probability  distribution  n(a)  of  radii  such  that 

n(a)da  ■  fraction  of  particles  with 
radii  in  (a,a+da) 


and 


amax 

min 


n(a)da 


1 


then  it  can  be  shown  that  the  volume  scattering  and  absorption 
coefficients  for  this  aerosol  are 


N  f 

aCr  J*-: 


max 

min 


oscaM  n (a) da 


(5.26) 


and 


“v.M  *  Naer  f™*  n«<ia 

amin 


(5.27) 


respectively.  The  phase  function  will  be 


Pv,M  ' 


fa  max 

Anin 


(i  +  i  )  n(a)da 
1  2 


hi.  “ 

min 


(ij  +  iz)  n(a)da 
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where  the  integral  over  dfl  *  dysd$s  is  an  integral  over  all 
scattering  angles;  its  appearance  in  the  denominator  guaran¬ 
tees  the  proper  normalization  of  P  ,  u  to  4ir  .  From  the  ex- 

V|M 

pression  for  o  in  terms  of  i  and  i  , 
r  sea  i  2 


°sca  "  T ?/,  [4 


and  from  Eq.  (5.26),  the  phase  function  may  be  written 


P"*  [i.tM.)  *  i.tK,)]  nCa)da  .  (5.28) 

V*M  */amin 


PR  and  Pv  M  may  be  combined  as  follows  to  yield  the  complete 
phase  function  for  scattering: 


0v,R  PR  +  Bv,M  Pv,M 


(5.29) 


where 


Bv,P.  +  Bv,M 


(5.30) 


is  the  total  volume  scattering  coefficient. 

Details  as  to  the  actual  computation  of  i^  ,  i^  , 

6  ,  and  P  are  omitted  here  for  brevity.  The  Mie  pro- 

V  ,M  V  ,M 

gram  has  been  coded  and  debugged  using  existing  tables  of  Mie 
functions,  and  every  effort  has  been  expended  to  keep  its  cost 
minimal,  in  view  of  the  fact  that  it  is  only  one  small  part  o* 
the  radiation  program.  An  elaboration  of  the  numerical  tech¬ 
niques  used  will  be  furnished  in  the  final  report. 
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5.2  BOUNDARY  CONDITIONS 

f  18) 

As  is  clear  from,  for  example,  Figure  1 . 1  of  Goody,'-  J 
the  spectrum  of  solar  radiation  and  the  spectrum  of  terrestrial 
radiation  overlap  hardly  at  all.  Therefore,  we  will  speak  of 
the  atmospheric  radiation  problem  a*,  two  separate  problems,  and 
discuss  the  boundary  conditions  for  each  problem  separately. 

5.2.1  Solar  Spectrum 

The  following  data  completely  specify  the  boundary 
conditions  for  the  solar  radiation  problem: 

•  solar  zenith  angle 

•  solar  constant  (flux  of  solar  energy  at  the 
top  of  the  atmosphere) 

•  solar  spectral  energy  distribution 

•  albedo  or  reflection  coefficient  at  the  sur¬ 
face  of  the  Earth. 

The  solar  zenith  angle  at  any  particular  location  on 
the  Earth’s  surface  is  a  function  of  the  day  of  the  year  and 
the  time  of  day.  Since  time  zones  are  so  irregular  as  to  be 
virtually  meaningless,  all  times  will  be  taken  as  Greenwich 
mean  times.  Then  the  solar  zenith  angle  0  er0°,90°l  will 

be  computed  from  the  four  variables  latitude,  longitude, 
calendar  date,  and  Greenwich  mean  time.  Leap  years  will  be 
accounted  for.  By  having  the  capability  to  specify  the  solar 
zenith  angle  in  this  fashion,  comparisons  can  be  made  between 
the  model's  predictions  and  experimental  data  gathered  at  any 
time  in  the  past. 

The  solar  constant  is  still,  surprisingly,  a  subject 
of  debate.  The  best  values  of  both  the  solar  constant  and  the 
solar  spectral  energy  distribution  seem  to  be  those  of 
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Thekaekara  and  his  collaborators^*^  which  will  be  used  in  the 
present  code.  The  solar  constant  will  be  adjusted  according 
to  calendar  date  because  of  the  varying  Earth-sun  distance, 
which  can  alter  the  solar  constant  percent  from  its  mean 
value.  Taking  the  origin  of  vertical  coordinates  z  *  0 
at  the  surface  of  the  Earth,  the  solar  boundary  condition 
can  be  expressed  mathematically  as 

■  sv  «(S  -  J!sun) 


for 


-l  <  y 


<  o 


0  <  <J>  <  2  it 


where  zq  is  the  vertical  coordinate  of  the  "top’'  of  the 
atmosphere  and  is  the  solar  intensity  at  frequency  v  . 

In  terms  of  the  azimuthally-averaged  intensity,  this  is 


Wu)  = 


Sv 

7? 


6(p  - 


ysun^ 


-1  <  V  <  0 


where 


sun 


=  -cos  0. 


sun 


The  sun  is,  of  course,  not  a  6-function  but  has  a  finite 
angular  width,  of  about  4°.  A  more  realistic  boundary  condi¬ 
tion  was  used  in  the  heating  rate  equation  for  no  scattering 
to  estimate  the  error  of  using  the  6-function.  The  conclusion 
was  that  the  absolute  error  produced  in  the  heating  rate  is 
always  inconsequential  (whereas  the  relative  error  in  the 
heating  rate  could  become  substantial  near  sunrise  or  sunset, 

v  >  0)  . 
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The  albedo  of  a  surface  is  the  ratio  of  the  outgoing 
to  the  incoming  flux.  Since  intensities,  not  fluxes,  are 
being  calculated,  the  specification  of  the  albedo  alone  is 
not  sufficient.  The  distribution  of  the  outgoing  or  reflected 
intensity  in  angle  is  also  required.  The  assumption  is  often 
made  that  the  Earth's  surface  is  a  "Lambert  reflector,"  mean¬ 
ing  that  the  reflected  intensity  is  isotropic  and  unpolarized 
regardless  of  the  angular  distribution  and  polarization  of  the 
incident  radiation.  Rough,  irregular  surfaces  approximate 
Lambert  reflectors.  If  the  incident  flux,  computed  from  the 
incident  intensities,  is  Fv  ^nc  ,  and  the  albedo  (which  may 
be  frequency-dependent)  is  Av  ,  then  for  a  Lambert  or  diffuse 
reflector 


C 

v,ref 


ylv(0,y)dy 


irlv(0,y) 


Av  Fv,inc 


Thus,  the  reflected  intensities  lv(0,y)  ,  0  <_  y  <_  1  ,  are 
specified  in  terms  of  the  albedo  and  an  integral  Fv  inc  of 
the  incident  intensities: 


IvC0,y) 


ylv(0,y)dy 


0  <  y  <  1 


Measurements  of  albedo  indicate  a  more  complicated 
situation  than  that  described  above.  On  cloudless  days, 
the  albedo  seems  to  be  fairly  constant  for  solar  zenith  angles 
®sun  —  »  an<*  t0  increase  rather  markedly  as  0gun  increases 

from  60°  to  90°.  In  particular,  this  phenomenon  is  observed 
for  the  ocean  and  for  desert  and  semi-desert  areas.  There¬ 

fore,  let  us  define  a  directional  albedo,  or  directional- 
hemispherical  reflectivity,  A^(y)  ,  which  is  the  reflected 
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flux  divided  by  (and  due  to)  an  incident  flux  irom  the  direc¬ 
tion  (0,4>).  (0e[O°,9O°]  is  the  angle  to  the  surface  normal.) 
Then  if  the  incident  intensity  from  direction  ft  is 
lv(0»n,<f>)  ,  '1  i  M  1  0  ,  the  incident  flux  will  be 

|li|  lv(0,p,$)  dQ 

which  will  cause  a  reflected  flux 

av C I u I )  |v|  iv(o,p,*)  dn  . 

Summing  over  all  incident  directions  leads  to  the  total  re¬ 
flected  flux 


p 

v  ,ref 


M  Av C  I W I )  IvC0,m  ,<*>) 


av(|m|)  Iv(0,m)  dp 


For  a  diffuse  reflector,  this  implies  a  reflected  intensity 
of 


fO 

Iv(0,M)  =2  J  |p|  Av(|y|)  lv(0,p)  dp 


0  <  p  <  1 


(5.31) 


Clearly,  the  directional  albedo  A^(p)  contains  no 

information  as  to  the  angular  distribution  of  the  reflected 

radiation.  Such  information  is  furnished  in  complete  detail 

(211 

by  the  bidirectional  reflectivity  p  , v  1  which  is  equal  to 
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the  reflected  intensity  Iv  ref(® »er »^r)  at  angles  0r  , 

4>r  ,  due  to  an  incident  intensity  Iy  inc(e><J0  at  angles 
0  ,  4>  ,  divided  by  the  flux  of  that  incident  intensity: 

„  fe  *  0  4  1.  "VrefC8’*’8.-’^ 

PvC9»*.8r.*r)  iVjjnc(6,*)  cose  an  • 

The  various  angles  are  defined  in  Figure  5.2.  If  the  reflec¬ 
ting  surface  is  isotropic,  as  the  Earth's  surface  largely  is, 
Pv  will  depend  only  on  the  difference  ~  4>r  , 

P.  “  P..(6  •  Since  I, _ r  is  of  differential  order 

v  v  r  r  v  ,rex 

with  respect  to  Iv  ^nc  (except  in  the  case  of  specular  reflec¬ 
tion)  ,  the  dfi  in  the  denominator  keeps  p  from  being  of 
differential  order.  The  factor  it  is  introduced  so  that,  if 
the  reflection  is  diffuse  (Iv  re£  independent  of  0r,  <J>r) ,  pv 
reduces  to  the  directional  albedo  Av  defined  above. 


Figure  5,2  -Geometry  of  reflection. 


The  total  reflected  intensity  at  angles  0r  ,  <J>r  is 
found  by  summing  Iv  rg£  over  all  possible  angles  0  ,  <J> 
of  incidence ; 
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VrefCW  ‘  Z  I»,r.£t°-*-,r-*P 


incident 

angles 


,  r2v  r-n/2 

j  J  J0  d6  sin0  pv^8 ,er’^"^r^ 


*  cose  . 


v,inc 


Phrasing  this  in  our  usual  notation, 


-  f2n  /*0 

Iv(0,y,4>)  -  £  d*’  J  i  dy'  |y'  |  Pv(|n'  I  »y,*'- 


♦) 


•  lv(0,y' ,*')  for  0  <  y  <  1  . 


Because  pv  must  be  periodic  in  its  third  argument, 


pv(y,y'  ,<|>+2ir)  =  Pv(y,y'  ,*) 


the  azimuthally-averaged  bidirectional  reflectivity 


_ ..  i  f  2* 

pv  ■  7?  J0  Kv(v,v',*-r) 


d<(» 


can  be  reduced  to 


>viy.y')  ■  ZJ  J  Pv(P,y’,J)  dj 


(5.32) 
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which  is  independent  of  .  This  result  makes  it  possible 
to  azimuthally  average  Eq.  (S.32)  to  yield 


lv(0,y) 


lu'l  P^dP'I.P)  lv(0,w*)  , 


0  <  p  <  1 


(5.33) 


This  is  the  most  general  reflective  boundary  condition. 

For  water  surfaces,  which  cover  about  three-fourths  of 

the  Earth,  computations  of  pv  are  possible  in  terms  of  the 

Fresnel  formulas  for  reflection,  the  index  of  refraction 

m  *  n  -in  of  water,  and  a  statistical  distribution  of 
1  2 

surface  slopes  (as  a  function  of  wind  speed).  The  work  of 
Chow^22^  is  exemplary  in  this  regard,  although  he  ignores  the 
imaginary  part  n^  of  the  index  of  refraction  in  the  IR  and 
uses  a  frequency-averaged  value  of  n^  in  his  computations. 

He  also  uses  the  simplest  analytical  approximation  to  Cox  and 
Munk's^23^  experimentally-determined  sea-slope  distributions. 
We  have  replaced  these  approximations  with  more  accurate  ones 
and  computed  tables  of  for  use  in  Eq.  (5,33).  The  tabu¬ 

lation  is  somewhat  simplified  by  the  reciprocity  relation  for 


Pv(e,er,<j>-<j>r)  =  Pv(er,e,<i>r-<f>) 


which  implies 


PV(P,P')  =  pv(p'  ,w)  • 

It  should  be  noted  that  only  surface  reflection  is 
accounted  for  in  the  above  computation  of  p^  .  Backscattering 


118 


3SR-795 


from  turbidity  (primarily  micro-organisms)  beneath  the  surface 

is  not  accounted  for,  although  measurements  in  the  Russian 

literature  indicate  that  this  effect  is  only  important  at  low 

(’241 

solar  elevations . v  1 

The  situation  vis-6-vis  reflectivity  data  for  land 
surfaces  is  much  worse  than  for  the  sea  surface  (cf.  Kondrat'yev, 
Ref.  24).  In  general,  only  albedos  are  available,  and  often 
not  even  as  a  function  of  frequency.  Therefore,  the  code  will 
have  several  options.  All  options  will  use  Eq.  (5.33);  how¬ 
ever,  if  only  directional  albedos  A^(p)  are  available,  dif¬ 
fuse  reflection  will  be  assumed  so  that  Eq.  (5.33)  reduces  to 
Eq ,  (5.31).  If  only  albedos  A^  are  available,  it  will  be 
presumed  that  Av(y)  ■  Av  ,  independent  of  y  .  And,  if  only 
frequency-averaged  albedos  A  are  available,  the  code  will 
take  A^  =  A  .  Thus,  as  improved  albedo  data  become  available, 
they  need  only  be  entered  into  the  appropriate  tables  and 
certain  option  flags  re-set. 

At  this  point  it  is  convenient  to  introduce  an  additive 
splitting  of  the  downward-directed  intensity: 

Tv  =  Tflar  ♦  if  «  ,  -1  <  y  <  0  .  (5.34) 

■fSolar  is  the  solar  beam  intensity,  and  is  the  "diffuse" 

intensity  produced  by  scattering  and  reflection  (thermal  emis¬ 
sion  being  practically  negligible  in  the  solar  spectrum).  The 
reason  for  introducing  this  splitting  is  that  the  solar  part 
of  the  intensity  is  essentially  a  6-function  in  angle  and  so 
exceedingly  difficult  to  represent  numerically.  The  remaining 
part  of  the  intensity,  T^1  ,  is  usually  smoothly  behaved  as 

a  function  of  angle  and  so  will  not  require  nearly  as  fine  an 
angular  mesh  for  its  representation  as  the  original  intensity 

I  would  have . 
v 
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The  solar  beam  intensity  can  be  found  from  the  trans¬ 
port  equation  in  which  only  extinction  processes  are  included: 


..ysolar 

^ - «vc*)i;olar 


1  <  p  <  0 


where  .  Imposing  the  boundary  condition 


I?.01" 


sun" 


1  <  p  <  0 


leads  to  the  solution 


ysolar 

Av 


<»*[£  £°  •'vCz')  dl’ 


-i  <  p  <  o  . 


Introducing  the  splitting  (5.34),  with  the  above  solution  for 
ysolar  ^  int0  the  full  transport  equation,  (5.9),  one  finds 

3ydiff 

p!j* _ «.(B  -  r11^) 

M3z  v'-v  J 


+  K 


\  f  °  Pv(z.u.n')  rjiff  U,V’)  dv 


(5.35) 


//  PvCz,u,u’) 


,  .  Tdif f 

dp'  -  !v 


e  s 

*  Pv^z,u,Msun^  exP 


r—  r°  Kv(z,)  dz' 

Hsun  J z 


for  -1  <  p  <  0  ,  and 
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I 


31 


VTT  "  ai(Bv  '  V 


*  6..I T  I  Pv(x,M,P')  T^iff  (X.u')  dll' 


*fo  •* 


1  _ 


Iv  (X.M'O  dll'  -  I, 


CS. 36) 


*  TT  exp 


-Lf\ 

^sun  •/ z  v 


(z*)  dz' 


for  0  _<  y  _<  1  . 

The  boundary  condition  on  at  the  top  of  the 

atmosphere  is  now  simply  a  homogeneous  one, 

-1  <  M  <  0  . 


rjiff  cx0,u)  •  o 


(5.37) 


The  reflective  boundary  condition  at  the  surface,  Eq.  (5.33), 
becomes 

Tv(0,u)  -  2  J  1m'  1  pv(|l,|u,|)  I^iff  (0 »M 1 )  dp' 


*  “?  I  “sun  I  Pv^'^sunU  exp 


1 

f  0  iev(x')  dz' 

M  i 

Hsun  J 

f0 

0  <  \i  <  1 


(5.38) 
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The  extra  terms  in  Eqs.  (5.35)  and  (5.36)  are  interpretable 
physically  as  scattering  sources  due  to  the  solar  beam;  the 
extra  term  in  Eq.  (5.38)  is  attributable  to  reflection  of  the 
solar  beam. 

A  final  boundary  quantity,  one  of  paramount  interest, 
is  the  solar  radiative  flux  into  the  ground,  which  determines 
the  heating.  It  is 


Fv(0)  -  2  n 


M)  dy 


ro 

2n  J  ^  yTJ1**  (0,y)  dy  +  2n 


f1  ^V0’ 

Jo 


y)  dy 


(5 


*  “sun  Su  exP 


-1 —  f 

ysun  Jo 


•cv(z')  dz’ 


When  this  is  integrated  over  frequency  and  added  to  the  total 
terrestrial  radiative  flux  out  of  the  ground,  the  resultant 
flux  determines  a  boundary  condition  for  a  ground  heating 
calculation.  If  the  flux  is  entirely  absorbed  within  the 
first  millimeter  or  so,  as  is  the  case  for  most  land  sur¬ 
faces,  then  it  determines  a  surface  heat  source;  if  it  pene¬ 
trates,  as  in  the  oceans  and  ice,  a  distributed  heat  source 
is  determined.  A  code  which  solves  the  heat  condition  equa¬ 
tion  with  prescribed  sources  has  been  developed  at  S3  during 
this  contract  period,  and  will  be  coupled  to  the  radiation 
code  for  studies  of  the  dynamic  interaction  of  the  radiation 
field  and  the  surface.  It  seems  likely  that  the  development 
of  this  coupled  radiation-surface  code  will  be  only  partially 
completed  by  the  end  of  the  present  contract  period,  and  will 
require  further  work  in  the  follow-on  period. 


.39) 
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5.2.2  Terrestrial  Radiation  Spectrum 

The  relevant  data  required  to  specify  the  terrestrial 
radiation  boundary  condition  are: 

o  the  surface  temperature,  T^ 

o  the  surface  emissivity  e^,  possibly  as  a 

function  of  the  angle  (to  the  surface  normal) 

8  of  emission. 

The  surface  temperature  T^  is  presumed  either  given 
as  an  input  variable  or  calculated  by  the  ground  heating  code 
discussed  in  Section  5.2.1. 

The  directional  emissivity  is  defined  as  the  ratio 
of  the  thermally  emitted  intensity  I v (Q  , ♦)  in  a  particular 
direction  to  the  black  body  intensity: 

ev(0,*)  ■  jVr  .  (5.40) 

VV  gJ 

We  shall  only  consider  isotropic  surfaces,  so  that  the  depend¬ 
ence  on  azimuthal  angle  <f>  disappears,  ev  =  ev(0)  or  ev  *  ev(y)  . 

If  the  emitted  intensity  Iv  is  independent  of  0  ,  ev  re¬ 
duces  to  the  more  familiar  hemispherical  emissivity  (ratio  of 
emitted  flux  to  black  body  flux  ttBv CT^)  ).  The  hemispherical 
and  directional  emissivities  have  been  measured  for  many 
kinds  of  surfaces.  The  angular  behavior  of  ev  is  similar 

for  all  electrical  non-conductors  (in  particular,  the  Earth's 
surface)  ;  it  is  nearly  constant  and  close  to  one  for  0  between 
0°  and  60° ,  then  it  falls  off  to  zero  as  0  increases  from  60° 
to  90°.  This  effect  has  never,  to  the  author's  knowledge, 
been  included  in  an  atmospheric  IR  radiation  model.  Since  the 
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IR  radiation  is  nearly  isotropic  (except  in  the  8-12  y  window 
region),  substantial  portions  of  it  approach  the  ground  at 
angles  of  60°  to  90°,  and  much  more  of  this  radiation  is  re¬ 
flected  than  a  constant  cv  model  would  indicate.  The  present 
model  will  incorporate  a  typical  angular  dependence  e(y)  for 
non-conductor  emissivities  and  use 

-  e<|0>  o(v)  (5.41) 

where  is  the  hemispherical  emissivity  for  the  surface 

in  question. 

The  directional-hemispherical  reflectivity  Av(y)  , 
the  bidirectional  reflectivity  Pv (y ,yr ,$-$r)  ,  and  its  azi¬ 
muthal  average  F..(y»y,)  were  defined  in  Section  5.2.1..  They 

v  r(211 
are  related  as  follows :v  J 


Av(y) 


yr  Pv(y,yr, ♦-♦r) 


dyr  yr  pv^y,yr^  * 


(5.42) 


Kirchoff's  Law  allows  us  to  relate  A  and  e  , 

v  v  * 

Av(y)  +  ev(p)  *  1  (5.43) 

f  211 

which  holds  without  restrictions.''  *  Hence,  if  the  reflec¬ 
tion  is  diffuse  so  that  Eq.  (5.31)  applies  for  the  reflected 
radiation,  the  surface  boundary  condition  may  be  formulated 
entirely  in  terms  of  e^(y)  » 
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IV(0.H) 


«v00  W  +  2 


•  [l  -  £v(l‘,'l>]  d>1' 


0  <  p  <  1  .  (5.44) 

The  first  (emission)  term  comes  from  the  definition  (5.40)  of 
emissivity . 

Only  for  water  surfaces  is  the  function  pv  obtain¬ 
able  in  the  IR.  It  is  theoretically  calculable  as  discussed 
in  Section  5.2.1.  From  it  we  may  obtain  A^(p)  ,  and  hence 
(p)  ,  according  to  Eqs.  (5.42)  and  (5.43).  For  such  sur¬ 
faces,  we  will  replace  the  second  term  in  Eq.  (5.44)  by  the 
exact  result  Eq.  (5.33),  eliminating  the  assumption  of  diffuse 
reflection. 

The  boundary  condition  on  the  terrestrial  radiation 
at  the  top  cf  the  atmosphere  is  homogeneous: 

Tv(z0,p)  *  o  ,  -1  _<  p  <  0  .  (5.45) 

5.3  DISCRETIZATION  OF  THE  TRANSPORT  EQUATION  FOR  NUMERICAL 
SOLUTION 

The  intensity  I  is  a  function  of  the  three  inde¬ 
pendent  variables  v  ,  z  ,  and  p  .  The  discretization  of  each 
of  these  variables  in  turn  is  discussed  in  Sections  5.3.1 
through  5.3.3.  The  actual  numerical  solution  of  the  transport 
equation  is  treated  in  Section  5.3.4. 

5.3.1  v- Discretization 

In  the  regions  of  the  spectrum  in  which  line  absorp¬ 
tion  is  important  (which  includes  most  of  the  terrestrial 
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radiation  spectrum  except  for  the  8-12p  "window"  region  and 
the  solar  spectrum  below  0.3p  and  above  lp)  the  absorption 
coefficient  varies  extremely  rapidly  with  frequency. 

So,  therefore,  will  the  intensity,  making  it  infeasible  to 
solve  for  Tv  at  a  set  of  discrete  v's  because  of  the 
large  number  of  v's  that  would  need  to  be  taken. 

Instead,  we  shall  solve  for  the  frequency- averaged 
intensities : 


IjU.lO  -  5 - i— -  fVl 

1  vi+l  vi  Ai,  v 


dv 


over  an  appropriately  small  number  of  frequency  intervals 
(vi»vi+P  *  ®ecause  this,  the  integro-differential  form 
(5.9)  of  the  transport  equation  is  unsuitable.  (It  is  not 
known  how  to  approximate  the  integral 


/.Vl  1  Tv  dV  • 


i 

in  which  both  and  Jv  are  violently  oscillating,  in 

terms  of  1^.)  We  must  use  the  integral  form  (see  Eq.  (5.2)  ) 


Iv(z,w)  -  lv(0,y)  exp 


\C  expf- 1 C,  dz' 


(5.46) 

dz" 


for  p  >_  0  ,  and 
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i  rz 

Iv(z,M)  =  Iv(zQ,y)  exp  -  -  J  kv(z')  dz' 


-  f  Jv(z",y)  exp  -  i  J  Kv(z»)  dz' 


dz"  (5, 


for  y  <  0  .  The  source  function  is 


Jv(z,y)  -  ct^(z)  Bv(T(z))  +  Sv(z,y) 


where 


A  Bv(z)  Z*1  —  — 

sv(z,y)  *  — 2 —  /  pvCz.y»v')  ivCz,n’)  dy' 


for  the  terrestrial  spectrum,  and 


sv(z,u)  ■  —2 


f°  PvU,u,w’)  Tjlff  (z,w')  dn' 


f1  - 

J  Pv(z,u,y')  Iv(z,y') 


dy 


S  _  l  -  rz 

■=—  Pfz,y,y  )  expl-  ~ —  /  0  k  C  z  * )  dz'! 

271  vv  ,M,Hsun'  (PSUn  j z  v  \ 


for  the  solar  spectrum.  (For  the  solar  spectrum,  I  should 

_i :  rr  v 

also  be  replaced  by  l“  in  Eq.  (5.47)  .)  Applying  the  fre 
quency  averaging  operator 


47) 
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to  both  sides  of  Eq.  (5.46)  , 

Ii(z,y)  =  ^(O.y)  Ti(y,0,z)  2^  (y,0,z) 


3T.  (y,z",z) 

(V,z”,z)  ypi -  dz" 


S^z’-.y)  Si  (y,z«’,z) 


(y,z”,z)  dz" 


(5.48) 


where  by  definition, 


(y  ,z  »z, 
r  12 


(5.49) 


Vu, 


i  +  1 


/Vl  exp 

'  v  J  2  °i(z)dz 

J*i  J 

dv 


(5.50) 


and 


a  B,(z)  rl 

Si(z,y)  *  —  J  Pi(z,y,y’)  I± (z ,y * )  dy’ 


(5.51) 
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for  the  terrestrial  spectrum,  while 


'  Bi(z)  /*0  diff 

>!<*•»>  ■  -V-  Jml  V*.*.*')  if 


(z , p  * )  dp1 


f  P.(z,p,p')  I. (z,p')  dp' 
Jn  1  1 


(5 


+  I?  Pi(z»P»PSun)  Zi  Cysun’zo'z)  Ti(ysun’zo'z) 


for  the  solar  spectrum.  S ^  is  taken  as  the  frequency- 
averaged  solar  flux 


s. - I -  P-1  s 

1  vi*l  -  vi 


dv 


since  this  is  the  form  in  which  solar  spectral  data  is  always 
presented.  The  quantities  3^  ,  ,  P^  are  the  correspond¬ 
ing  quantities  8V  ,  ,  P^  evaluated  at  the  mid-point 

of  the  frequency  interval, 


.52) 
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vt  -  >i(vi  +  v.+1)  . 

The  of  Eq.  (5.50)  are  called  transmission  functions. 

An  important  approximation  has  had  to  be  made  in  de¬ 
riving  Eq.  (5.48),  aside  from  the  relatively  trivial  one  of 
approximating  the  slowly  varying  (in  frequency)  functions  , 

Bv  ,  and  by  their  values  at  the  mid-point.  It  is  the 
commuting  of  the  frequency-average  with  lv(0,y)  ,  and  with 
Tv(z,y')  in  the  scattering  term,  and  the  replacement  of  these 
quantities  by  1^(0, y)  and  I^(z,y')  .  In  view  of  boundary 
conditions  such  as  Eq.  (5.44),  the  boundary  intensities 
Ty(0,y)  may  be  rapidly  varying  functions  of  v  (unless  ev  =  1)  . 

Similarly,  in  the  presence  of  substantial  line  absorption, 
Iv(z*m')  will  vary  rapidly  with  v  .  This  hurdle  has  re¬ 
sulted  in  two  divergent  bodies  of  independent  research  in 
atmospheric  radiation;  one  school  neglects  scattering  (e.g., 

Kyle,  Ref. 26  ) t  the  other  neglects  absorption  (e.g.,  Sekera, 

Ref.  27).  The  primary  thrust  of  the  first  school  is  the 
accurate  calculation  of  the  transmission  functions  T.  .  Once 

A  1 

these  are  known,  and  ■  0  ,  »  1  ,  the  numerical 

integration  of  Eq.  (5.48)  is  almost  trivial.  The  second  school 
concerns  itself  with  techniques  for  solving  the  integral 
equation  (5.48)  when  T^  ■  1  ,  ■  0  (in  which  case  the  com¬ 

mutation  of  the  frequency  average  with  lv(0,y)  and 
Iv(z,y')  is  a  valid  approximation).  To  the  author's  know¬ 
ledge,  no  one  in  all  the  vast  literature  on  atmospheric  radia¬ 
tion  has  considered  simultaneously  line  absorption  and  scatter¬ 
ing. 

The  reason  for  this  apparent  lacuna  is  that,  over 
large  portions  of  the  solar  and  terrestrial  spectrums ,  either 
scattering  or  line  absorption  is  dominant.  They  might  only 
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by  comparable  in  magnitude  in  the  near  infrared  region  (l-5p) 
where  there  are  some  weak  bands,  and  in  parts  of  the  8-12y 
window  region.  Therefore,  the  error  that  we  make  in  doing  the 
frequency-average  of  the  scattering  term  in  Eq.  (5.46)  will 
tend  to  be  large  in  only  a  small  fraction  of  the  frequency 
intervals;  presumably  these  errors  will  have  little  impact 
on  the  frequency- integrated  flux  of  Eq.  (5.14). 

A  potentially  serious  approximation  is  the  replace¬ 
ment  of  1^(0, y)  by  1^(0, y)  ,  particularly  in  the  strong  IR 
absorption  bands.  If,  in  Eq.  (5.44),  the  emissivity 
indeed  falls  to  zero  over  a  span  of  angles  of  30°  or  so,  then 
the  second  term  in  that  equation  is  not  negligible.  Since 
the  incident  intensities  1^(0, y')  in  that  term  will  be 
rapidly  varying  functions  of  v  ,  so  will  the  reflected  in¬ 
tensities,  and  hence  1^(0, y)  itself.  A  mitigating  circum¬ 
stance  in  favor  of  the  approximation  is  that,  in  the  strong 
IR  bands,  the  surface  boundary  condition  will  become  unim¬ 
portant  after  about  the  first  kilometer  or  two;  that  is,  the 
transmission  function  T^(y,0,z)  multiplying  1^(0, y)  be¬ 
comes  negligibly  small  for  z  >  1-2  km.  Nevertheless,  one 
would  hope  to  do  the  boundary  layer  correctly.  Therefore, 
ways  to  skirt  this  difficulty  are  actively  being  sought. 

One  possible  avenue  of  approach  is  to  revert  to  the 
integro-differential  form  (5.9)  and  attempt  to  deal  with  the 
integral 


I  dv 


(the  f requency -averages  of  the  other  terms  are  trivial). 
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Defining 


fitl  “v 

Jv. 


\  dv 


P*1  rv 

•'v . 


dv 


the  integral  may  be  approximated 


(5.53) 


1  rVi 
Ev  J 


IV  dv 


1  -  e 


-hvi/kT 


)“i 


cu  is  a  function  of  both  z  and  \i  ,  in  general.  It  also 
depends  on  the  intensity,  of  course,  which  is  the  source  of 
the  difficulty.  Nevertheless,  by  calculating  for  vari¬ 

ous  typical  intensity  fields  in  the  atmosphere  (obtained  by 
detailed  calculations)  regularities  might  emerge  which  would 
allow  us  to  pick  a  universal  a^(y,z)  .  If  this  were  possible, 
it  would  not  only  alleviate  the  difficulties  discussed  above 
but  would  actually  be  simpler  to  tabulate  than  ,  which 

depends  on  three  arguments. 

The  computation  of  the  transmission  functions  T^  hos 
a  long  history.  The  earliest  attempts  were  based  on  band 
models,  in  which  simple  analytical  representations  of  line 
strengths,  positions,  and  shapes  were  assumed.  As  accurate 
line-by-line  absorption  data  has  become  available ,  both 
from  theory  and  experiment,  transmission  function  computations 
have  incorporated  it.  Such  detailed  line-by-line  transmission 
function  computations  are  incredibly  expensive  in  terms  of 
computer  time.  Considering  that  sometimes  integration  steps 
as  small  as  0.001  cm‘‘  must  be  taken^^  ,  and  that  the  region 
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of  significant  absorption  extends  from  10,000  cm"1  (ly)  to 
250  cm"1  (40y)  ,  with  only  a  few  gaps,  the  magnitude  of  the 
problem  becomes  apparent.  As  an  example,  Kyle^2^  used  15 
minutes  of  CDC  6600  time  to  compute  transmission  functions 
between  a  single  pair  of  atmospheric  levels  z^  and  z^  for 
a  single  value  of  y  ,  and  for  the  wavelength  interval 
1.7y  -  20y.  Multiply  this  by  the  number  of  angles  and 

by  the  number  of  pairs  of  levels  JjN  (N  -1)  ,  and  the  computing 
time  to  obtain  a  complete  set  of  transmission  functions  be* 
comes  truly  formidable  (for  6  angles  and  15  levels  it  would 
be  1 5 7*s  hours).  Fortunately,  transmission  functions  are  not 
terribly  sensitive  to  the  temperature  profile*,  and  so  could 
be  tabulated  once  and  for  all  for  various  standard  profiles 
(tropical,  mid -latitude,  polar,  for  summer  and  winter,  for 
example).  This  would  restrict  us,  however,  to  always  using 
the  same  pressure  levels  and  angles,  which  could  be  a  large 
disadvantage . 

In  view  of  the  expense,  in  terms  of  both  human  and 

computer  time,  of  generating  transmission  functions  from 

scratch,  we  have  decided  to  take  advantage  of  the  scheme  of 

McClatchey,  et.al.,^30^  for  generating  these  functions. 

It  falls  into  the  category  of  an  empirical  fit  to  known  data, 

and  is  the  most  sophisticated  in  a  long  line  of  such  empirical 

fits.^3*^32^  McClatchey  has  used  detailed  line-by-line 

absorption  data  to  compute  transmission  functions  for  20  cm*1 

f  k) 

intervals,  then  has  fit  them  with  empirical  functions  fv  J 
such  that 


* 

Kyle,  private  communication. 
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The  superscript  k  refers  to  molecular  species;  there  are 
separate  f's  for  the  uniformly  mixed  absorbers  (CO2 ,  N20, 
CH4,  CO,  02,  N2) ,  for  water  vapor,  and  for  ozone.  The  single 
(kl 

argument  Aw^'  is  an  attempt  to  sum  up  the  information  con¬ 
tained  in  p,z  ,  and  z  into  a  single  "effective  absorber 
1  2 

amount"  along  the  slant  path  in  question.  It  is  calculated 
according  to  empirical  prescriptions  which  best  fit  the  real 
data;  the  variable  mixing  ratios  of  H20  and  0^  are  taken 
into  consideration  in  these  prescriptions.  The  total  trans¬ 
mission  function  is  the  product  of  the  individual  ones 


which  is  an  approximation  also,  but  an  excellent  one  according 
to  several  authors  .  ^ The  functions  f^)  are  tabulated 
in  a  subroutine  called  LOWTRAN  which  we  have  obtained  from 
McClatchey  and  implemented  on  our  computer. 

While  McClatchey's  transmission  functions  will  be 
used,  it  would  be  desirable  in  the  longer  range  to  develop  a 
code  which  could  generate  its  own  transmission  functions  di¬ 
rectly  from  the  raw  absorption  data.  Among  the  advantages 
of  this  are: 


(1)  frequency  intervals  Av  could  be  chosen 
as  desired. 

(2)  improvements  in  the  raw  absorption  data 
could  be  readily  incorporated;  and 

(3)  the  approximation  of  an  "effective  absorber 
amount"  could  be  replaced  by  a  more  accurate 
one  (the  Curtis -Godson  approximation^^  for 
example) . 
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We  have  already  had  discussions  with  both  Kyle  and  McClatchey 
about  the  possibility  of  obtaining  their  absorption  data  and 
transmission  function  generators.  They  have  both  mentioned 
the  high  cost  in  terms  of  computer  time  to  generate  trans¬ 
mission  functions  from  scratch,  but  in  the  interests  of  a 
definitive  radiation  calculation  we  believe  these  costs 
would  be  justified.  However,  for  the  present  time  we  have 
deferred  these  discussions  in  order  to  pursue  other  aspects 
of  the  code  development. 

5.3.2  y-Discretization 

The  angular  variable  y  appears  both  as  a  parameter 
and  as  a  variable  of  integration  in  the  transport  equation 
(5.9)  or  Eqs .  (5.46)  and  (5.47).  Because  measurements  and 
theoretical  computations  all  indicate  that  the  terrestrial 
radiation  intensity  and  the  diffuse  part  of  the  solar  inten- 
sity  (I?1  t)  are  fairly  smooth  functions  of  angle,  these 
intensities  may  be  represented  by  their  values  at  a  rela¬ 
tively  few  angles  y .  ,  i  ■  1 ,  . . . ,  N  .  In  order  to  do  the 
flux  integrals,  Eq.  (5.11),  the  intensity  will  be  assumed  to 
vary  in  a  piecewise-linear  fashion  (I  ■  1^  ♦  yl^^)  between 
the  points  yi  at  which  it  is  calculated.  To  ensure  consis¬ 
tency,  the  scattering  source  integrals 

Si(Zj,yk)  =  sijk  *  dU' 


will  be  done  under  the  same  assumption, 


N  -1 

y 


E 


V1 

E 


/Vl  *  “''ij.n)  d"' 

yn 


o(0)  .CO)  „  Cl)  j(i) 

qijk,n  i , j ,n  qijk,n  i,j,n 
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where  the  moments 


of  the  phase  function  are  defined  as 


q 


(a) 

ijk,n 


Ua  pi(zj  »Wk»P)  dy  . 


The  moments  q^  ,  q^  are  to  be  computed  and  stored  be¬ 
fore  the  calculation  of  the  intensities  begins.  It  is  perhaps 
worth  noting  that  will  have  to  be  computed  over  a  finer 

angular  mesh  than  y^  in  order  to  ensure  accuracy  in  the 
numerical  integrations  leading  to  q^  ,  q^  . 

An  exception  to  the  statement  that  varies 

smoothly  in  angle  occurs  when  there  is  substantial  aerosol 
scattering.  Aerosol  scattering  notoriously  produces  a  strong 
forward  peak  in  the  scattered  intensity.  This  forward  peak 
is  almost  as  troublesome  numerically  as  the  solar  beam  it¬ 
self  (which  we  eliminated  by  the  splitting  of  Eq.  (5.34)  )  be¬ 
cause  it  necessitates  a  dense  mesh  of  angles  around  the  solar 

angle  y„..„  .  Therefore,  we  shall  use  a  method,  tested  by 

sun 

Hansen, in  which  the  radiation  scattered  into  a  narrow 
forward  cone,  say  ±2°  about  the  forward  direction,  is  re¬ 
garded  as  unscattered.  Mathematically,  this  amounts  to 
truncating  the  forward  peak  from  the  phase  function  P.,(z»wQ) 
and  decreasing  the  scattering  coefficient  0v  accordingly. 

For  dust  and  haze,  which  are  optically  thin,  this  seems  like 
a  reasonable  procedure;  in  clouds,  on  the  other  hand,  there 
might  be  a  cumulative  error  after  many  scatterings  which  is 
not  small.  Hansen  shows  this  not  to  be  the  case,  however, 
and  so  the  method  is  fully  viable  for  all  types  of  aerosols. 


5.3.3  z-Discretization 

Pressure  coordinates  Pj  will  be  substituted  for 
elevations  z^  in  the  model  by  making  the  hydrostatic 
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assumption, 


dp  *  -pg  dz 

This  will  facilitate  comparison  with  the  GCM's  and  agrees 
with  general  meteorological  practice. 

There  is  no  discretization  of  the  pressure  coordinate 
which  is  ideal  in  all  regions  of  the  spectrum.  The  fundamen¬ 
tal  criterion  for  vertical  zoning  is  that  the  source  function 
Jv  may  not  change  substantially  from  zone  to  zone.  If  it 
did,  interpolations  necessary  to  do  the  Jv  integrals  in 
Eqs .  (5.46)  and  (5.47)  would  be  too  inaccurate.  In  the 
infrared,  this  means  that  BV(T)  (and  therefore  T)  and  the 
transmission  function  T^  may  not  change  substantially  be*> 
tween  any  two  zone  centers.  Since  6°/km  is  a  typical  lapse 
rate,  zones  should  probably  not  exceed  2  km  in  width.  Kyle, 
in  hJs  IR  model^^,  used  15  2-km  levels  surmounted  by  a 
sixteenth  level  from  30  km  to  the  top  of  the  atmosphere.  In 
most  of  the  solar  spectrum,  on  the  other  hand,  a  much  coarser 
vertical  resolution  could  be  tolerated,  since,  except  within 
aerosol  layers,  the  scattering  source  varies  considerably 
less  rapidly  with  height  than  BV(T)  , 

'  hese  considerations  suggest  that  a  dynamic  assign¬ 
ment  of  zoning  structure  would  be  highly  desirable,  based  on 
an  examination  both  of  the  source  function  and  the  spectral 
interval  involved.  In  the  case  of  the  scattering  source 
function,  which  involves  the  intensity,  we  would  examine 
the  scattering  coefficient  8V  and  phase  function  Py  in¬ 
stead.  In  absorption-dominated  spectral  intervals,  we  will 
zone  so  that  the  relative  change  in  temperature  T  and  trans¬ 
mission  function  T^  from  zone  to  zone  will  be  bounded.  In 
scattering-dominated  spectral  intervals,  we  will  zone  so  that 
the  relative  changes  in  Bv  and  Pv  (at  selected  pairs  of 
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angles  y  ,  y')  are  similarly  bounded.  In  spectral  intervals 
where  absorption  is  comparable  to  scattering,  both  T  ,  T\ 
and  0  ,  P  will  be  required  to  be  bounded  in  their  zone-to- 

zone  variations. 

Aerosol  layers,  particularly  clouds,  will  be  much 
better  resolved  by  such  a  zoning  scheme  than  they  would  be  by 
problem-independent  schemes  (e.g.,  fixed  2-km  levels).  By  the 
same  token,  zones  will  not  be  wasted  in  regions  across  which 
very  little  is  happening  to  the  radiation  field. 

The  problem  with  such  a  dynamic  zoning  scheme  is  that 
the  fluxes  F^  in  different  frequency  intervals  i  will  not 
be  calculated  at  the  same  levels.  By  the  very  nature  of  the 
dynamic  zoning  scheme,  however,  the  individual  fluxes  F^Cz^) 
will  not  vary  greatly  between  z^  and  zj+i  »  therefore,  a 
polynomial  interpolation  scheme  (a  parabolic  fit,  for  example) 
can  furnish  F^(z)  at  any  intermediate  level  between  z.  and 
zj+l  •  Hence,  the  various  fluxes  F^Cz^)  can  be  reduced  to 
common  levels  z^  prior  to  summation  over  frequency. 


FUk) 


This  dynamic  zoning  scheme  is  compatible  with 
McClatchey's  transmission  function  method  discussed  in  Sec¬ 
tion  5.3.1.  It  would  also  be  compatible  with  exact  trans¬ 


mission  functions,  precomputed  for  specific  levels  z 


CO) 


for  such  functions  would  vary  smoothly  enough  with  level 
that  they  could  be  interpolated  to  the  levels  of  interest. 


5,3,4  Numerical  Solution 

The  basic  method  for  the  numerical  solution  of  the 
transport  equation  will  be  discrete  ordinates  with  scattering 
iteration . 
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The  discrete  ordinates  method  was  first  proposed  by 
Chandrasekhar.^^  It  involves  fixing  the  angle,  y  =  , 

in  the  integral  form  (5.46)  and  (5.47)  of  the  transport  equa 
tion  and  integrating  from  one  boundary  to  the  other  in  dis¬ 
crete  steps  Az  .  The  method  is  simply  ray-tracing,  account 
ing  for  sources  and  sinks  along  the  ray  trajectory  (and, 
incidentally ,  ignoring  the  slight  bending  of  the  ray  due  to 
refraction) . 

The  Planck  function  Bv(T(z))  in  Eqs .  (5,46)  and 
(5.47)  is  to  be  interpolated  linearly  in  z-space  between 
points  z.  at  which  it  is  known,  for  purposes  of  doing  the 
z-integratior.s  in  (5.46)  and  (5.47)  numerically.  This  en¬ 
sures  that  the  diffusion  limit  of  radiative  transport  will 
be  recovered.  To  see  this,  suppose  there  is  no  scattering, 
that  =  a ^  is  independent  of  z  ,  and  that  we  are  suffi¬ 
ciently  many  optical  mean  free  paths  from  the  boundary  that 
the  boundary  term  in  Eq.  (5.46)  is  negligible.  Then 


Iv(z»P) 


U 


rz  /  \ 

o* 

/„ 

1  exp 

-  jp  (2-2") 

•  . 

dz' 


A  partial  integration  leads  to 


I v  Cz ,  m) 


bw(tcu)  -  bv(t(<»)  exp 


-  -a'  z 
V  v 


JC 


z  3B 


TPT  exp 


$*i(*-*") 


dz' 


We  neglect  the  second  term  on  the  right-hand  side  because 
the  optical  path  from  the  boundary  a^z  is  assumed  large. 
In  the  third  term,  if  B^  is  linear  in  z"  ,  then  3Bv/3z" 
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is  constant,  and 


Bv(T(z))  - 


w 

7 


3B 

v 

TT 


This  is  the  diffusion  approximation. 

The  scattering  iteration  process  can  best  be  ex¬ 
plained  by  writing  Eq.  (5.46)  schematically  as 


*  L<V 


Into 

involving 
on  I  .  and 


IQ  we  have  lumped  all  the  known  terms  (the  ones  not 


V- 


Miv) 


is  a  double  integral  operator  operating 
is  the  scattering  source  term.  The  iter¬ 


ation  we  shall  use  to  solve  this  equation  is  si 


7(n*l) 

v 


I0  ♦  LCI^0) 


where  the  initial  iterate  is 


t(0) 


( 37) 

This  process  is  known  to  always  converge. v  '  Furthermore, 
when  there  is  a  significant  amount  of  absorption  the  con¬ 
vergence  is  rapid.  On  the  other  hand,  if  absorption  is 
negligible,  as  for  visible  radiation  in  a  water  cloud,  and 
there  are  many  mean  free  paths  of  scattering,  convergence  is 
painfully  slow.  This  can,  to  some  extent,  be  circumvented 
by  initializing  the  intensity  in  the  scattering- thick  region 
not  from  IQ  ,  but  from  one  of  the  two  famous  analytic  approxi 
mations  which  are  attached  to  the  names  of  Eddington  and 
Schwarzschild,  respectively.^®^  A  variant  of  this  technique 
has,  in  fact,  been  successfully  applied  at  S*.^9^ 
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The  radiative  transfer  code  will  be  included  in  a 
realistic  boundary  layer  code  that  has  been  developed  under 
this  contract.  This  boundary  layer  code  is  described  in  Ap 
pendix  E.  The  coupled  code  will  provide  an  opportunity  to 
evaluate  radiative  effects  in  comparison  to  ordinary  thermo 
dynamic  effects  (latent  heat  transfer,  etc.)  in  the  atmos¬ 
phere  . 
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6.  FUTURE  STUDIES 

This  section  of  the  report  will  attempt  to  outline 
further  code  development,  modifications,  and  application  of 
these  codes  to  test  problems  which  will  be  attempted  during 
the  next  six  months  of  this  research  contract. 

6.1  HAIFA  CODE  DEVELOPMENT  AND  APPLICATIONS 

Investigations  now  in  progress  will  be  completed. 
These  include: 

(1)  Further  test  and  modification  of  the  ver¬ 
sion  of  HAIFA  incorporating  compressibility, 

(2)  Further  test  of  the  version  of  HAIFA  in¬ 
corporating  moisture  effects. 

(3)  Further  modification  of  the  Poisson  solver 
to  include  non-cyclic  inlet  and  outlet 
boundary  conditions. 

Several  modifications  to  the  HAIFA  codes  will  be  made 
in  order  to  study  phenomena  which  arc  not  yet  understood  com¬ 
pletely.  Among  the  major  code  developments  will  be: 

(1)  the  addition  of  a  turbulence  scheme  to 
the  basic  HAIFA  code,  and 

(2)  the  addition  of  Coriolis  terms  to 
determine  their  importance  in  mountain 
wave  problems  on  the  meso-scale. 
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It  is  also  desirable  to  take  account  of  effects  in  three 
spatial  dimensions.  It  will  not  be  possible,  however,  to 
develop  a  code  and  carry  out  such  calculations  under  the 
current  contract,  the  scope  and  magnitude  of  the  problem 
will  be  investigated  for  future  consideration. 

Additional  calculations  to  be  performed  during  the 
remaining  period  of  the  current  contract  include: 

(1)  Investigation  of  the  effect  of  mountain 
shape  on  lee  waves ,  and 

(2)  Calculation  of  a  problem  using  a  fully 
compressible  code  for  comparison  with 

a  standard  problem  calculated  with  HAIFA. 

These  results  and  those  from  previously  calculated  cases  will 
be  analyzed  more  quantitatively  to  characterize  the  momentum 
and  energy  transports.  As  a  part  of  this  analysis,  the 
perturbed  pressure  field  will  be  calculated  through  solution 
of  the  governing  Poisson  equation. 

The  final  two  months  of  the  contract  period  will  be 
spent  in  analyzing  and  attempting  to  parameterize  the  results 
of  the  many  calculations  in  such  a  manner  as  to  be  useful  in 
the  global  circulation  model  calculations. 
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APPENDIX  A 

DERIVATION  OF  BOUSSINESQ  EQUATIONS 

The  conservation  equations  governing  macroscopic  fluid 
motion  are  frequently  simplified  for  problems  of  thermal  con¬ 
vection  by  introducing  certain  approximations  which  are  at¬ 
tributed  to  Boussinesq.  These  approximations  can  best  be 
summarized  by 

(1)  fluctuations  in  density  which  appear 
with  the  advent  of  motion  result  prin¬ 
cipally  from  thermal  (as  opposed  to 
pressure)  effects,  and 

(2)  in  the  conservation  equations  of  mass 
and  momentum,  density  variations  may 
be  neglected  except  when  they  arc 
coupled  to  the  gravitational  accelera¬ 
tion  in  the  buoyancy  force. 

These  approximations  are  examined  in  the  derivation  of  equa¬ 
tions  presented  below. 

The  general  equations  of  mass  and  momentum  conserva¬ 
tion  are 


-pV*V 


(A.l) 
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pdV  _  .yp+  y.^  .  pgJ  .  (A.  2) 

For  purposes  of  this  derivation  the  viscous  stress  tensor  P 
will  be  dropped  from  the  equations.  The  equation  of  state 
will  be  assumed  to  be  of  the  form 


P  ■  P(P»T)  . 


(A. 3) 


The  basic  approximation  to  be  made  may  be  examined  by  the  fol¬ 
lowing  procedure: 

(1)  Let  f  represent  any  one  of  the  state  variables. 

It  will  be  expressed  in  the  following  form 

f  *  fm  ♦  VZ>  +  f CAM) 

where 

fB  «  space  average  of  f 

f  (Z)  ■  variation  of  f  in  the  absence 
of  motion 

f'(x,z,t)  ■  fluctuations  in  f  resulting 

from  fluid  motions. 

(2)  If  a  scale  height  is  introduced  as 


H(f) 


1  dfo 

IT 

m 


-1 


(A.  5) 


the  basic  approximation  is  that  the  fluid  be  confined  to  a 
layer  whose  thickness,  d  ,  is  much  less  than  that  of  the 
scale  height  (d  <<  H)  . 
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In  particular,  Eq.  (A. 5)  implies  that  d/H(p)  <<  1.  On  inte¬ 
grating  this  latter  condition  over  the  layer,  one  concludes 
that 


Ap 

S— -  5  e  «  1  ,  (A. 6) 

Hm 

where  ApQ  is  the  maximum  variation  of  pQ  across  the  layer. 

It  is  also  required  in  non-linear  investigations  to 
make  the  additional  restriction  that  the  motion  induced  fluc¬ 
tuations  do  not  exceed,  in  order  of  magnitude,  the  static 
variation,  i ,e . , 


£l 

pm 


<  0(e) 


(A.  7) 


Condition  A. 7  must  be  verified  a  posteriori  from  solutions 
of  the  problem.  In  the  absence  of  motion  and  introducing 
Eq.  (A. 4),  the  vertical  component  of  Eq.  (A. 2)  is 


^o 

5T~  "  ’gpm  '  gpo 


(A.  8) 


Introducing  the  hydrostatic  relation  into  Eq.  (A. 2),  we  have 

P(|y  +  v*Vv)  *  -Vp’  -  gp'  k  ,  (A. 9) 

We  may  introduce  Eqs.  (A. 4)  and  (A. 6)  into  the  continuity 
Eq.  (A.l)  to  obtain 


V.»v 
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Hence  to  order  e  ,  Eqs .  (A. 9)  and  (A. 10)  may  be  written 


||  *  V-V7  - 


-  r  v 

Hm 


p'  C 
k 

Ko 


(A. 11) 


V«v  =  0  (A. 12) 

A 

In  Eq.  (A.  11)  we  have  retained  the  term  ge(p*/Ap0)  k  even 
though  it  contains  e  as  a  factor.  This  procedure  is 
necessary  if  we  are  to  study  convection  problems  in  the 
Boussinesq  approximation,  and  the  following  justification  may 

be  made:  The  quantity  measures  the  characteristic 

acceleration  of  the  fluid.  Now  the  system  is  driven  by 
fluctuations  of  the  density  field,  and  hence  we  must  insist 
that  the  characteristic  acceleration  be  of  order  (gcp'/Ap^  . 
This,  in  turn,  forces  the  conclusion  that  the  acceleration 
of  gravity  is  always  much  greater  than  the  characteristic 
acceleration,  i.e., 
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APPENDIX  B 

WATER  PRODUCTION  TERM 

The  water  productioi  term  consists  of  three  physical 
phenomena  which  can  add  to,  subtract  from,  or  change  the  state 
of  the  water  in  the  atmosphere.  Thus  the  production  term  Pr 
is  written  as  the  sum  of  (.1)  the  evaporation  of  rainwater  out¬ 
side  of  the  cloud,  0(r-rs)  ;  (2)  the  conversion  of  cloud 
water  to  rainwater;  and  (3)  the  accretion  of  cloud  water  by 

rainwater,  S  ; 

'  a 


Pr  -  BCr-rsl  ♦  Sc  +  Sa 

Cb 

where 

3 

is  the  evaporation  parameter  and  is 
assumed  constant, 

Sc 

is  a  linear  function  of  the  cloud 
water  content;  and 

S 

a 

is  a  variable  dependent  on  both  the 
cloud  water  content  and  the  rainwater 
content , 

f  71 

The  expressions  for  SQ  and  Sa  ,  taken  from  Orville v  ' 

are : 

S  *  aft  -  i  ) 
c  L  c  crJ 

CB.2) 

Sa  ■  4.6  *  10'3  lc(lr)0,95  . 

CB.3) 
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Orville  has  run  test  problems  using  ranges  of  the  constant 

a  from  10"*  sec'*  to  2  *  10'^  sec*  ana  values  of  1 _  from 

-4  -1  cr 

0  to  5  x  10  gn  gin 

Variations  of  these  parameters  and  their  values  will 
also  be  studied  as  part  of  the  research  on  determining  the 
moisture  effects  on  mountain  lee  waves  and/or  the  drag  force 
on  the  air  flow. 
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APPENDIX  C 
EDIT  QUANTITIES 

The  momentum  flux  (wave  drag)  associated  with  gravity 

waves  is  fundamentally  different  from  other  known  momentum 

transport  processes  like  surface  frictional  drag  inasmuch  as 

fCll 

it  may  act  across  deep  atmospheric  layers.  Sawyerv  * 
pears  to  be  the  first  to  have  pointed  out  that  in  stratified 
flow  the  pressure  is  systematically  higher  on  the  upstream 
side,  resulting  in  a  drag  force  on  the  obstacle,  and  a  cor* 
responding  drag  of  opposite  sign  on  the  air  stream. 

For  steady  flow  over  an  obstacle,  Vergeiner ^2) 
shows  that  the  equation  for  horizontal  momentum,  Eq.  (2.1) 
may  be  integrated  over  a  slab  between  the  mountain  and  a  fixed 
height  H  to  give  a  drag  on  the  mountain  equivalent  to 


D(drag) 


(cuk)z„h  d*  - 


I 


H 

(P  ♦  PU2) 

obstacle 

height 


x-L 

dz 

x*  -L 


(C.l) 


where  L  and  -L  indicate  lengths  upstream  and  downstream 
from  the  obstacle.  If  (p  ♦  pu2)  is  the  same  for  upstream 
and  downstream  (L  •*  »)  ,  the  momentum  flux  puw  dx  is 
constant  as  a  function  of  z  and  equal  to  the  drag.  In  the 
linearized  case,  the  drag  may  be  transformed  by  using  the 
linearized  equation  for  horizontal  momentum  into 
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-p(z) 


dx 


u ' w '  dx 
o  o 


(C .  2) 


where  the  o  subscripts  indicates  values  taken  at  the  surface 
of  the  mountain. 

The  problem  we  are  concerned  with  in  the  research  is 
the  calculation  of  a  drag  for  the  air  flow  over  a  mountain 
where  the  numerical  calculations  are  not  capable  of  being  run 
to  steady  state  values.  The  quantities  presently  being  edited 
from  the  results  of  the  numerical  calculations  are 


2L  I  u'w'  dx  or  u*w' 

J-L 

where  2L  is  the  distance  of  the  horizontal  grid.  The  quanti 
ties  u’w'  have  been  obtained  by  two  different  methods. 
Initially,  an  averaging  method  was  used,  i.e.,  u'  was  found 
from 


u' 


(ui.j 


♦  u. 


)/2  -  u 


(C.3) 


where 


I 

“  •  r£  <uij  *  Vi,j)/2  • 

i-1 

v*  was  found  as  (v.  .  ♦  v.  .  .)/ 2  .  In  this  manner, 
the  cell-centered  velocities  were  found  and  then 
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1 

2U 


dx 


was  calculated  by  averaging  the  product  of  the  two  perturba¬ 
tions  velocities  over  the  horizontal  grid  length. 

The  second  method  of  finding  the  drag  consisted  of 
editing  Crowley's  second  order  scheme  to  find  the  vertical 
flux  of  horizontal  momentum.  Crowley's  scheme  for  the  momen¬ 
tum  flux  is  written  in  finite  difference  form  is 


1 

7 


r  .  .  - 

i  ,j*l 


[(ui.j-i  *  uij)- 


At 

II 


(C.  4) 


By  averaging  this  quantity  over  the  horizontal  grid  the  result 
is  just  equal  to  the  product  of  the  horizontal  and  vertical 
perturbation  quantities.  This  follows  by  noting  that 


I 


i 


♦ 


.  (C .  5) 


A  comparison  of  these  two  methods  shows  essentially 
identical  results. 

The  results  as  described  above  for  the  problems  com¬ 
pleted  to  date  are  given  in  Section  4  of  this  report.  In  order 
to  attempt  to  understand  these  results,  various  lengths  of  the 
horizontal  average  were  used  in  obtaining  values  of  u 'w'  . 
These  results  are  also  discussed  in  Section  4. 
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APPENDIX  D 
HAIFA  CODE  LISTING 

A  listing  of  the  HAIFA  code  which  is  presently  oper¬ 
ating  on  the  UNIVAC  1108  at  Systems,  Science  and  Software  (S1) 
is  included  in  this  appendix.  The  advection  scheme  is  second 
order  only  in  this  version  of  the  code.  Working  versions  at 
S’  allow  either  second  order  or  fourth  order  schemes. 
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APPENDIX  E 

A  COMPUTER  CODE  FOR  THE  ONE -DIMENSIONAL 
BOUNDARY  LAYER 

E.l  INTRODUCTION 

The  radiative  transfer  task  of  the  Climatology  con¬ 
tract  calls  for  the  development  of  an  accurate  numerical 
scheme  and  application  cf  it  to  the  thermodynamics  of  the 
atmosphere  and  soils.  Ultimately,  this  scheme  is  to  be  used 
to  perform  calibration  calculations  of  the  radiative  sub¬ 
routines  of  General  Circulation  Models  of  the  Earth-’s  Climate. 

In  order  to  develop  a  realistic  radiative  transfer 
code  system  it  is  also  necessary  to  take  account  of  effects 
which  strongly  influence  the  thermodynamic  state  of  the  at¬ 
mosphere.  By  virtue  of  its  strong  influence  on  long  wave 
radiation  it  is  important  to  take  account  of  water  vapor  and 
cloud  moisture  in  the  atmosphere.  In  addition,  the  state  of 
the  lower  atmosphere  is  strongly  influenced  by  turbulent  trans¬ 
fer.  Consequently,  we  have  formulated  and  are  testing  a  1-D 
computer  code  to  evaluate  changes  in  the  atmosphere  resulting 
from  the  effects  of  radiative  transfer,  Coriolis  force,  turbu¬ 
lent  momentum,  heat  and  moisture  transfer,  and  subsidence. 

E.2  FORMULATION 

The  formulation  depends  only  on  the  vertical  coordi¬ 
nate  and  corresponds  to  an  atmosphere  in  which  all  properties 
are  horizontally  homogeneous,  The  atmosphere  is  described  by 


El 
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the  eastward  horizontal  velocity  component  u  ,  the  northward 
velocity  component  v  ,  the  temperature  T  ,  and  the  relative 
humidity  q  ,  each  of  which  may  depend  on  the  vertical  coordi¬ 
nate  z  .  The  atmosphere  is  assumed  to  be  instantaneously 
in  hydrostatic  equilibrium  but  changes  in  temperature  result 
in  a  vertical  subsidence  velocity  w  .  Pressure  gradients 
in  the  horizontal  directions  are  also  taken  into  account  but 
they  are  assumed  to  depend  only  on  z  , 

The  1-D  description  of  the  atmosphere  and  the  hydro¬ 
static  approximation  permit  the  use  of  a  Lagrangian  formula¬ 
tion  in  which  the  atmospheric  pressure  is  a  convenient  measure 
of  the  mass  of  the  atmosphere  above  the  mass  element  in  ques- 
tion.  The  altitude  above  the  surface  z  and  the  pressure  p 
are  related  by 


where  the  density  p  is  to  be  determined  from  the  equation 
of  state,  and  pQ  is  the  atmospheric  pressure  corresponding 
to  the  surf are  z  ■  0  , 

The  equation  for  the  transient  boundary  layer  have 
been  treated  by  many  authors,  e.g.,  Estoque,^*^  Pandolfo,^^ 
and  Sasamori.^^  We  choose  to  employ  the  temperature  as 
dependent  variable,  rather  than  potential  temperature,  because 
temperature  is  more  closely  related  to  the  radiative  properties 
which  form  the  most  important  aspect  of  this  investigation. 
Using  the  pressure  independent  variable  the  equations  are: 


E2 


&  -  *(*-*,>  ♦  (vg)  • 

3T  “-f (u-ug)  *  g'pjp  (vl?)  . 

-  (If  *  «}r  +  «*p|p  (ktp||)  -  *Pr^X  , 

Quantities  appearing  in  Eq.  (E.2)  are  given  by: 


QUANTITY 

SYMBOL 

EQUATION 

REMARKS 

density 

P 

■  It 

determined  by  equation 
of  state 

vertical 

velocity 

w 

-  ir 

3t 

Lagrangian  derivative 
of  altitude 

Coriolis 

parameter 

f 

■  2  flsin$ 

<p  is  latitude;  0  is 
angular  velocity  of 

Earth 

Geostrophic 

wind 

Ug 

>  .  1  |£ 
Tp  3y 

y  is  north-south  distance 

Vg 

"  IP  If 

x  is  east-west  distance 

Adiabatic 
lapse  rate 

r 

g  is  gravitational  con¬ 
stant;  C  is  specific  heat 
of  air  ”  at  constant 
pressure 

Radiation 

flux 

F 

to  be  determined  in  radia¬ 
tive  subroutine 

Turbulent 

diffusivity 

K 

to  be  determined  in 
k-subroutine 

E3 
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The  indicated  time  derivatives  m  Eq.  (E.2)  are  to  be  formed 
at  constant  pressure;  consequently,  they  are  Lagrangian  deri¬ 
vatives  in  that  they  evaluate  changes  associated  with  a  par¬ 
ticular  air  mass  element.  The  advection  term  associated  with 
vertical  motion  is  included  in  these  terms. 

E.3  EDDY  DIFFUS1VITY 

The  coefficients  of  turbulent  transfer  which  appear 
in  Eq.  (E.2)  has  been  developed  through  a  combination  of 
theoretical  considerations  and  empirical  observations,  A 
number  of  expressions  for  these  quantities  are  available 
representing  different  weightings  of  the  data  and  greater  or 
lesser  sophistication  in  incorporating  theoretical  considera¬ 
tions. 

Several  of  the  investigators  assume  that  the  four  eddy 
coefficients  are  equal  in  the  two  momentum  equations  and  the 
temperature  and  relative  humidity  equation.  Such  is  the  case 
in  the  work  of  Sasamori^E3^  who  uses  the  eddy  diffusion  coeffi¬ 
cients  developed  by  Yamamoto  and  Shimanuki . ^E4^  The  same 
assumption  is  made  by  Estoque,  et.al.^E^  who  attribute  their 
expression  to  Blackadar . ^E5^  In  our  current  work  we  use  the 
prescriptions  of  Pandolfo^E2^  who  has  modified  the  Monin- 
Obukhov  formulae  as  presented  by  Kitaigorodsky . Pandolfo 
takes  account  of  differences  between  the  coefficients  for 
momentum  exchanges  and  those  of  heat  and  moisture.  He  also 
imposes  limitations  on  the  magnitude  of  the  coefficients  cor¬ 
responding  to  the  case  of  extreme  stability. 

The  expressions  for  the  exchange  coefficients  are; 

Inversion  Conditions  (Ri  >  0) 


Kv  "  *T  "  Kg  “  k2(Z+Z0)2  |t  Cl+oRi) 


(E.3) 
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Lapse  Forced-Convection  Conditions  (-0.048  <  Ri  £  0) 


KU  “  Kv  *  k’Cz+Zo>2  JZ  U-om* 


(E.4) 


KT  -  Kq  -  Ku  (1-aRi) *2 


Lapse  Free-Convection  Conditions  (Ri  <  -0.048) 

Kt  -  Kq  -  h(z.z0)>jf(|I  ♦  r)  I** 


(B.5) 


where 


K  -  K  -  K  l  '1/2  Ri!“1/6 
*u  *v  *T  c  R1| 


a  ■  -3,0 

c  -  3(|)^  (^)J  1 
h  -  (0.4)2  (|)3/Z 

Ri  ■  f[H  * r  *  °-61T  • 

The  above  values  of  K  are  to  be  restricted  in  range 
in  order  to  avoid  unrealistic  conditions.  According  to  Pandolfo, 
the  K-values  should  lie  within  the  following  ranges: 


104  <  K  <  107  cm2/sec 


z  >  100m 


102  <  K  <  107  cm2/sec 


z  <  100m  . 
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m  0 


ii 

•m  • 


T 


m  * 


Computed  values  falling  outside  of  the  above  range  are  replaced 

< 

by  the  adjacent  bounding  values. 

Clearly,  these  exchange  coefficients  are  incapable  of 
taking  into  Account  such  effects  as  penetrative  convection  and 
advection  of  turbulence,  being  formulated  on  the  assumption  of 
steady  state  conditions.  We  hope  to  be  able  to  take  account 
of  these  effects  in  the  future.  The  influence  of  penetrative 
convection  has  been  estimated  by  Deardorff^7^  and  Estoque.^E®^ 
Dynamic  effects  of  turbulence  have  been  treated  in  varying 
degrees  of  rigor  by  Pritchett  and  GawAin,^E^  Harlow,  et.al,,^1^ 
and  by  Donaldson. ^E**^ 

E.4  DIFFERENCE  EQUATIONS 

The  Eq.  (E.2)  are  solved  as  a  set  of  coupled  difference 
equations  in  time  and  space.  The  difference  formulation  must 
satisfy  requirements  of  accuracy,  stability,  and  computational 
efficiency.  Several  considerations  affecting  accuracy,  stability 
and  efficiency  are  discussed  below. 

Large  shear  of  the  horizontal  wind  in  the  lower  portion 
of  the  atmospheric  boundary  layer  results  in  rapid  transport 
by  turbulence  and  correspondingly  large  transient  adjustments 
in  response  to  perturbations  of  the  boundary  layer.  In  order 
to  take  account  of  this  turbulent  transport  in  a  computationally 
efficient  way  we  have  formulated  the  equations  implicitly;  the 
result  is  an  unconditionally  stable  numerical  integration 
scheme  which  permits  the  time  interval  to  be  chosen  in  accord 
with  accuracy  considerations.  The  alternative  explicit  formu¬ 
lation  imposes  the  requirement  that  the  time  interval  satisfy 
the  inequality 


At 


-nr(ip) 
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In  the  lower  boundary  layer,  the  time  interval  permitted  by 
the  above  expression  will  be  very  small;  the  diffusivity  K 
is  large  and  the  desired  pressure  interval  Ap  will  be  small. 
The  implicit  formulation  which  we  have  selected  (discussed 
below)  requires  some  additional  calculations  to  solve  the 
sets  of  simultaneous  linear  equations.  However,  the  result 
is  a  system  of  equations  which  aro  stable  for  very  large  time 
intervals.  Time  intervals  are  determined  almost  entirely  by 
considerations  of  accuracy. 

The  changes  in  the  wind,  temperature  and  relative 
humidity  are  concentrated  predominantly  in  the  lower  layers 
of  the  boundary  layer.  From  the  standpoint  of  the  accuracy 
and  efficiency  of  the  numerical  integration  it  is  very  de¬ 
sirable  to  introduce  more  zones  in  the  region  of  rapid  change 
near  the  ground  than  higher  in  the  atmosphere  where  much 
smaller  changes  occur.  In  order  to  achieve  spatially  vari¬ 
able  resolution  with  accuracy  it  is  necessary  to  consider 
carefully  the  difference  formulation.  In  the  following 
scheme  the  difference  equations  retain  second  order  accuracy 
in  regions  of  variable  spatial  intervals.  The  resulting 
system  of  equations  is  capable  of  representing  the  atmosphere 
accurately  through  the  use  of  a  finely  resolved  layer  in  the 
lower  boundary  layer  and  increasingly  coarse  resolution  in 
the  higher  atmosphere. 

We  now  consider  the  difference  equations  correspond¬ 
ing  to  Eqs .  (E.2).  The  equations  are  to  be  solved  for  the 
primary  dependent  variables,  u  ,  v  ,  T  ,  q  on  a  discrete 
mesh  on  both  of  the  indpendent  variables,  pressure  p  and 
time  t  .  The  resulting  difference  equations  are  to  be 
solved  by  marching  the  solution  forward  by  successive  incre¬ 
ments  At  of  the  time.  In  order  to  use  a  closed-form  solu¬ 
tion  of  the  implicit  system,  we  linearize  those  terms  of  the 
equations  which  depend  on  the  advanced  time. 
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We  denote  a  discrete  value  of  the  time  by  superscript 
n  ,  i.e.,  tn**  -  tn  +  Atn  .  The  atmospheric  pressure  is  par¬ 
titioned  into  intervals  Ap^  such  that  ♦  Ap^  ■  P^+Jj  • 

The  difference  equations  corresponding  to  Eqs.  (E.2)  are: 


u»+1  .  un 

^/..n+1  \  .  _n  ^n  i+1  ui 

‘(vi  -  \,i)  —  -  — 


„n+l  n+1 
u,  -  u.  , 
,n  i  i-1 


‘  (KuC)i-L  — - —  ,  (B.6) 

u  1  11  4Pi  ♦  APi.jJ 


vn+1  -  vn 
vi  vi 


(ug,i  ‘  ui+1)  *  °i  (Kvp)i+>j 


v**1  .  vn+1 

l+l  vi 

APi+i  ♦ 


n+i  _  n+1 

,»n  vi  vi-l 

'  ^  vP^ i-4  - 

*  Ap.  +  Api_1 


,  (H.7) 


Ti  *  li 


J*U  ■  fU  >\  °L/„n  „n  V 

\  Ap,  7  JIrrT,i*%  KT,i-»»y 

Tn+1  Tn+1  Tn+1  Tn+1 

_n  _,n  *i+l  ‘  li  rv  ^n 
+  o.  (K-rPJi+L  "  (•'•rP J  l 

H  Api+1  ♦  ApA  *  APi  ♦  APi_x 
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n+1  _n 

2l_ Li..; 

At  1 


_  _n+l  „n+l 

(y)J+%  :  qi 


Ap 

i*i  *  *p 

_n+l 

„n+l 

<»i 

‘  «i-l 

Apt 

*  4Pi-i. 

(E.9) 


In  the  above  equations  we  have  defined 


2g*P? 

TpT 


and  the  diffusion  coefficients  are  to  be  evaluated  from 
Eqs.  (E. 3)  through  (E.S)  using  appropriate  centered  difference 
representations  of  the  dependent  variables  from  cycle  p. 

The  density  is  obtained  from  the  equation  of  state  as 

follows : 


where 


.n+1 


RT 


l 

n+1 


Pi  *  ‘stPi-Hj  *  PiV  • 

The  altitude  corresponding  to  the  pressure  also  depends  on 
time  by  virtue  of  the  hydrostatic  readjustment  of  the  vertical 
column  under  the  changing  temperature: 
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where  I  is  the  maximum  value  of  i  corresponding  to  the 
zone  adjacent  to  the  ground. 

The  vertical  velocity  is  obtained  as  a  difference 


where 


.n+1 


zi  zi 


At 


n 


«*!♦*  *  *i-%> 


Considering  the  above  as  a  system  of  simultaneous  equations 
for  the  unknown  quantities  u”+1  ,  v”+1  ,  T^+1  ,  g"+1  ,  we 

note  that  the  equations  are  linear  and  are  uncoupled  in  the 
following  way:  The  T-equation  and  q-equation  are  not  coupled 
to  each  other  or  to  the  u  or  v  equations;  the  u  and  v  equations 
are  coupled  together  through  the  Coriolis  terms.  Consequently, 
the  T  and  q  equations  can  each  be  represented  as  a  tri¬ 
diagonal  equation  for  a  scalar  quantity.  The  u  and  v  equa¬ 
tions,  however,  are  conveniently  represented  together  as  a 
tridiagonal  system  of  equations  for  a  vector  quantity  having 
thfe  two  components  u^  and  v^  . 

The  scalar  equations  can  be  represented  in  the  form 

+  Ci  *i-l  "  Di  ‘ 


ElO 


Ai  *i*l  +  Vi 


(E.10) 
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For  the  T-equation, 


Tn+i 

1 


and 


gj  AtCKTp)j.ns 
*pin  +  ^ 


4Pi  ♦  APi.l 


Di 


r  At 


¥i+h  fi-*i 


_n 
o . 


APi 


'  wi  "  l|  (ICT,i+Js  '  1CT,i-Ji) 


+  T. 


For  the  q-equation, 


.  n+1 

*i  ■  Oi 


and 


Ai 


Api+1  ♦  Apj 
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Bi 


-  1 


Ci 


* 


oj  ^t(Kap)?.^ 

APj  ♦  APi_! 


The  vector  equation  can  also  be  represented  in  the  form  of 
Eq.  (E.l)  where 


is  a  vector  quantity  and  the  coefficients  are  matrices  having 
the  form 


> 
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Ci 


,Ui  -  Atf  V 


Di 


(  ‘  ‘'I 

' v”  ♦  Atf  u  . • 

1  g»i 


where 


U 


♦ ,  °i  Atcy)j^ 

APi  ♦  Ap.+1 


v+ .  °i 

4pitl  *  LVl 


,n 


U 


^(kup)?.% 

♦  4Pi_i 


°j  *t(Kyp)?.% 

4Pi  ♦  iPi-j 


and  I  is  the  identity  matrix. 

All  of  these  systems  of  equations  can  be  solved 
readily  by  Gaussian  elimination.  Taking  advantage  of  the  tri¬ 
diagonal  form  of  the  equations,  the  solution  algorithm  re¬ 
duces  to  evaluation  of  coefficients  recursively  in  one  forward 

and  one  backward  sweep  through  the  mesh.  The  algorithm  for 

fE12) 

vector  equations  is  discussed  by  Richtmyer  and  Morton. 


E.5  BOUNDARY  VALUES 

The  calculational  region  extends  from  the  ground,  where 
the  pressure  has  the  assumed  ground  level  hydrostatic  value,  to 
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an  arbitrary  altitude  having  a  specified  pressure.  Boundary 
conditions  are  required  at  the  top  and  bottom  of  the  mesh  to 
close  the  system  of  equations.  We  have  not  investigated  these 
conditions  carefully  (they  will  be  affected  further  by  the 
radiative  treatment  at  the  ground),  but  are  using  the  follow¬ 
ing  set:  At  the  ground  level  the  velocity  is  zero,  correspond¬ 
ing  to  the  viscous  boundary  condition,  the  temperature  has  a 
specified  value,  and  the  relative  humidity  is  given  the  satura¬ 
tion  value  corresponding  to  the  ground  temperature.  At  the 
top  of  the  mesh  the  velocity  takes  the  geostrophic  value  and 
temperature  and  humidity  are  specified. 
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